---
title: "Simulate HLM with Mixtures of Normal Random Effects"
author: "Carolyn J. Anderson"
date: "2/7/2023"
output:
  html_document:
    toc: yes
---
# Introduction

The R covered here reproduced mini-study on impact of distribution of $U_j$s on results.  This corresponds to the R script title which is given below.

The function that simulates empty/null HLM where distribution of $U_{oj}$
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.

As written is simulates 4 cases by varying $\sigma^2$: 1, 9, 25 and 100

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

You can change values of $\sigma^2$.

# Set Up 

Load lmer & lmerTest
```{r, libraries}
library(lme4)
library(lmerTest)
```

Set working directory
```{r, wkdictory, eval=FALSE}
setwd("D:\\Documents\\Dropbox\\edps587\\lectures\\7 inference2")
```

Read in data
```{r, hsbData}
hsb <- read.table("hsball.txt",header=TRUE)
```

A quick look at the data
```{r, hsbOK}
head(hsb)
```


# Function to simulate very simple HLM

The function is named 'simmix' and takes as input the following:

*   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 
parameter estimates of $\tau_{00}$.  

The definition of the function 

```{r, fundef}
simmix <- function(N, nj, s, tau, inmu) {
   k <- 1                                # index for micro
   Uoj <- matrix( ,nrow=N,ncol=1)        # initalize matrix for U_0j
   yij <- matrix( ,nrow=N*nj,ncol=1)     # initalize matrix for y_ij
   cut  <- runif(N,min=0,max=1)          # draw uniform cut-points
   for(j in 1:N){
     mu <- inmu                          # read in value of mu
     if (cut[j]<0.5) {                   # use cut-point
        mu <-  -inmu                     # if <0.5, then mu=-inmu
     }                                     # i.e., draw of other distribution
     Uoj[j,1] <- rnorm(1,mean=mu,sd=tau) # draw random U_oj
     id <- j                             # index for macro
     for (i in 1:nj){ 
       yij[k,1]  <- Uoj[j,1] + rnorm(1,mean=0,sd=s) # compute y_ij
       k <- k+1                          # up-date index for mirco
     }
   }
   id <- as.data.frame(rep(1:N,each=nj))   # finishing up

   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
}
```

# Examples 

All of these examples set $N=1000$ groups, $n_j=5$, $\tau_{00}=1$ and $\mu= \pm 2$.

First level of $\tau_{00}$
```{r}
N <- 1000
nj <- 5
t <- 1
mu <- 2
```

What is varied in this set is the value of the within standard deviation
```{r}
s <- 1
```
You can change any of the above, but you might also want to change the titles in the graphs below to correspond to your input values.



Simulate data, graph distribution of $Y_{ij}$ and $U_{0j}$and fit HLM:
```{r}
simmix(N=N ,nj=nj, s=s, tau=t, inmu=mu)
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)
```

Change $\sigma^2=1$ to $\sigma^2=9$,
```{r}
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)
```

Changing $\sigma^2=25$,
```{r}
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)
```

...and now $\sigma^2=100$,
```{r}
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)
```

Let look at histograms of all historgrams of estimated $\hat{U}_{ij}$ in a single figure.  Note that to get the $\hat{U}_{ij}$s we need to first get them by using coef(model.x), which produces a "mer" object of random effects
```{r}
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") ))
```


Another way to look at these distribution and how close they are to normal is to use qqplots.
```{r}
par(mfrow=c(2,2))

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") ) )
```