#
# R for lecture on random effects
# Fall 2017
# C.J.Anderson
#
#  Function that simulates empty/null HLM where distribution of Uoj
#    is a mixture of two normal distributions.  There are a number
#    of things that you can change easliy.
#
#  After the function is are some examples of it's use....
#
# Simulates 4 cases by varying sigma2: 1, 9, 25 and 100

#    o Plots Yij and Uoj from the simulated data
#    o Fit models to data
#    o Plots historgram of estimated uoj from models
#    o Plot qqplots of estimated uoj from models


# load lmer & lmerTest
library(lme4)
library(lmerTest)

# Set working directory
setwd("C:\\Users\\cja\\Documents\\Dropbox\\edps587\\lectures\\7 inference2")

# Read in data
hsb <- read.table("hsball.txt",header=TRUE)

# look oK?
head(hsb)


#########################################################################
# study on effect of dist U_oj in very simple model (empty/null HLM)
#########################################################################
# Function to simulate very simple HLM
#
# N:  number of macro units
# nj: number of micro per macro
# s:  sigma (standard deviation of within group
# tau: tau (standard devations between groups)
# inmu:  mean of the two mixtures for Uoj (i.e., +/- inmu)
#
# You can change the input to see how even a little violation messes up 
#  paramter estimates of tau^2.  
########################################################################

simmix <- function(N,nj,s,tau,inmu) {
   k <- 1
   Uoj <- matrix(0,nrow=N,ncol=1)
   yij <- matrix(0,nrow=N*nj,ncol=1)
   cut  <- runif(N,min=0,max=1)
   for(j in 1:N){
     mu <- inmu
     if (cut[j]<0.5) {
        mu <-  -inmu
     } 
     Uoj[j,1] <- rnorm(1,mean=mu,sd=tau)
     id <- j 
     for (i in 1:nj){
       yij[k,1]  <- Uoj[j,1] + rnorm(1,mean=0,sd=s)    
       k <- k+1
     }
   }
   id <- as.data.frame(rep(1:N,each=nj))

   simdata <- cbind(k,id,as.data.frame(yij))
   names(simdata) <- c("index", "id", "yij")
   simdata <<- simdata

   Uoj <- as.data.frame(Uoj)
   names(Uoj) <- c("Uoj")
   Uoj <<- Uoj
}

par(mfrow=c(2,2))

t <- 1



simmix(N=1000,nj=5,s=1,tau=t,inmu=2)
simdata1 <-simdata
hist(simdata1$yij,breaks=15,xlab=expression(Y_ij), main=expression(paste("model.1: ",sigma^2,"=1,",  tau, "=1, ICC=0.500") ) ) 
Uoj1 <-Uoj 
hist(Uoj1$Uoj,breaks=15,xlab=expression(U_0j), main=expression(paste("model.1: ",sigma^2,"=1,",  tau, "=1, ICC=0.500") )  )
model.1 <- lmer(yij ~ 1 +(1 | id), simdata1,REML=FALSE)
summary(model.1)



simmix(N=1000,nj=5,s=3,tau=t,inmu=2)
simdata2 <- simdata
hist(simdata2$yij,breaks=15,xlab=expression(y_ij),main=expression(paste("model.2: ",sigma^2,"=9,",  tau, "=1,   ICC=0.100") ))
Uoj2 <-Uoj
hist(Uoj2$Uoj,breaks=15,xlab=expression(U_0j),main=expression(paste("model.2: ",sigma^2,"=9,",  tau, "=1,   ICC=0.100") ))

model.2 <- lmer(yij ~ 1 +(1 | id), simdata2,REML=FALSE)
summary(model.2)


simmix(N=1000,nj=5,s=5,tau=t,inmu=2)
simdata3 <- simdata
hist(simdata3$yij,breaks=15,xlab=expression(y_ij),main=expression(paste("model.3: ",sigma^2,"=25,",  tau,"=1, ICC=0.0385") ) )
Uoj3 <-Uoj
hist(Uoj3$Uoj,breaks=15,xlab=expression(U_0j),main=expression(paste("model.3: ",sigma^2,"=25,",  tau,"=1, ICC=0.0385") ) )

model.3 <- lmer(yij ~ 1 +(1 | id), simdata3,REML=FALSE)
summary(model.3)

simmix(N=1000,nj=5,s=10,tau=t,inmu=2)
simdata4 <- simdata
hist(simdata4$yij,breaks=15,xlab=expression(y_ij),main=expression(paste("model.4: ",sigma^2,"=100,", tau,"=1, ICC=0.0099") ) )
Uoj4 <-Uoj
hist(Uoj4$Uoj,breaks=15,xlab=expression(U_0j),main=expression(paste("model.4: ",sigma^2,"=100,", tau,"=1, ICC=0.0099") ) )

model.4 <- lmer(yij ~ 1 +(1 | id), simdata4,REML=FALSE)
summary(model.4)

par(mfrow=c(2,2))


## For historgrams of 
#   coef(model.x) produces a "mer" object of random effects, but the following
#     put random effects into form that can be used to plot, etc.
#

par(mfrow=c(2,2))

random1.uoj <- coef(model.1)$id[,"(Intercept)"]
hist(random1.uoj, main=expression(paste("model.1: ",sigma^2,"=1,",  tau, "=1,   ICC=0.500") ) )

random2.uoj <- coef(model.2)$id[,"(Intercept)"]
hist(random2.uoj,main=expression(paste("model.2: ",sigma^2,"=9,",  tau, "=1,   ICC=0.100") ))

random3.uoj <- coef(model.3)$id[,"(Intercept)" ]
hist(random3.uoj,main=expression(paste("model.3: ",sigma^2,"=25,",  tau,"=1, ICC=0.0385") ))

random4.uoj<- coef(model.4)$id[,"(Intercept)" ]
hist(random4.uoj,main=expression(paste("model.4: ",sigma^2,"=100,", tau,"=1, ICC=0.0099") ))




# qq-plots of random effects
plot(ranef(model.1),main=expression(paste("model.1: ",sigma^2,"=1,",  tau, "=1, ICC=0.500") ) )
plot(ranef(model.2),main=expression(paste("model.2: ",sigma^2,"=9,",  tau, "=1, ICC=0.100") ) )
plot(ranef(model.3),main=expression(paste("model.3: ",sigma^2,"=25,",  tau,"=1, ICC=0.0385") ) )
plot(ranef(model.4),main=expression(paste("model.4: ",sigma^2,"=100,", tau,"=1, ICC=0.0099") ) )