# Edps 587
# Spring 2019
# c.j.anderson
#
# 1. Simulate data from an a random intercept model with 2 levels
# 2. Fit and compare models fit to random intercept data (too simple, 
       correct model and too complex.
# 3. Simulate data from a random intercept and slope model
# 4. Fit and comparge models fti to intercept and slope model
#
# These could be writen more efficiently but they work and
# can be relatively easy to modify
#
# Using screenreg from the texreg package can save you a lot
# of typing.  To get output from texreg into word (and keep 
# spacing, use
#    o "courier new" font in word
#    o landscape (especially if have lots of models).
#          portrait works better for a few models and lots 
#          effects in them
#    o single space with 0 pt between lines
#    o edit model titles to be more informative
#    o IF needed use smaller font (my least favorite thing to do).
#
library(lme4)
library(texreg)

#################################################
# 1. simulated random intercept data            #
#################################################

# Set parameters for simulation:
N <-  100         # number of clusters/groups = 100 * 50
nj <- 20          # number of observations for group j

# Set fixed effects parameters
gamma00 <- 5
gamma01 <- 2
gamma10 <- 3

# Set variance of random intercept, slope, covariance, and level 1 residual
tau2 = 2
sigma2 = 4

# This will hold the simulated random intercept data
R.intercept <- matrix(999,nrow=N*nj,ncol=7)

# Need keep track of row in simulated data
index <- 1

###  Now simulate the random intercept data 
# (Outer Loop) Sample the marco units
    for (macro in (1:N)) {
     U0j= sqrt(tau2) * rnorm(1) 
     age=50 * runif(1) + 20      # this makes age range between 20 - 70
	 
#   (Inner loop) Within a marco unit, take a sample of micro units;
       for (micro in (1:nj)) {

# Create level 1 vexplanatory variable (draw from standard normal)
		xij = rnorm(1) 

# Create level 1 residual with sigma2= 2	
		z = rnorm(1)   
		Rij=sqrt(sigma2)*z;  

# Computed response/outcome variable		
        yij = gamma00 + gamma01*xij + gamma01*age + U0j + Rij
		R.intercept[index,1:7] <- c(macro,micro,yij,age,xij,U0j,Rij)

		index <- index + 1
	 }  # <--- end inner loop, go to next micro unit; 

} # <--- end outer loop, go to next macro unit;

R.intercept <- as.data.frame(R.intercept)
names(R.intercept) <- c("macro","micro","yij","age","xij","U0j","Rij")

#####################################################
#  2. fit models to simulated random intercept data #
#####################################################
# Correct model
summary( model.1 <- lmer(yij ~ xij + age + (1 | macro),data=R.intercept, REML=FALSE) )

# Overfit model
summary( model.2 <- lmer(yij ~ xij + age + (1 + xij | macro),data=R.intercept, REML=FALSE) )

# Too simiple 
summary( model.0 <- lmer(yij ~ xij  + (1 | macro),data=R.intercept, REML=FALSE) )

# Compare with each other and with the parameters used to simulate
screenreg(list(model.0,model.1,model.2))



######################################################
# 3. Simulate random intercept and slope model       #
######################################################

# Set parameters for simulation:
N <-  100         # number of clusters/groups = 100 * 50
nj <- 20          # number of observations for group j

# Set fixed effects parameters
gamma00 <- 5
gamma01 <- 2
gamma10 <- 3

# Set variance of random intercept, slope, covariance, and level 1 residual
varU0j = 4
varU1j = 1
covU0jU1j= .4
sigma2 =  4

 
# To get desired covariance matrix of U0j and U1j, use cholskey root
T.mat <- matrix(c(varU0j, covU0jU1j, covU0jU1j, varU1j),nrow=2,ncol=2)

w <- chol(T.mat)

# This will hold the simulated values
R.intslope <- matrix(999,nrow=N*nj,ncol=8)

# Need because I'm not efficient
index <- 1

# Simulate the data
# (Outer Loop) Sample the marco units
    for (macro in (1:N)) {
	
      U0j <- w[1,1]*rnorm(1)
      U1j <- w[1,2]*U0j + w[2,2]*rnorm(1) 
	  
      zj = 50 * runif(1) + 20      

# (Inner loop) Within a marco unit, take a sample of micro units;
       for (micro in (1:nj)) { 
      
	    xij=rnorm(1)
	   
        z=rnorm(1);   
        Rij=sqrt(sigma2)*z;      

# Computed response/outcome variable		
        yij = gamma00 + gamma10*xij + gamma01*zj + U0j + + U1j*xij + Rij

	    R.intslope[index,1:8] <- c(macro,micro,yij,zj,xij,U0j,U1j,Rij)

		index <- index + 1
	 }  # <--- end inner loop, go to next micro unit; 

} # <--- end outer loop, go to next macro unit;

R.intslope <- as.data.frame(R.intslope)
names(R.intslope) <- c("macro","micro","yij","zj","xij","U0j","U1j","Rij")



# Correlation matrix of random effects that were simulated --- is this what we want?
x <- R.intslope[ , c(6,7,8)]
cor(x)

################################################
# 4.  Fit Models to simulated data             #
################################################

# Fit model too simple
summary( simple <- lmer(yij ~ xij + zj + (1  |macro), data=R.intslope, REML=FALSE) )

# Correct model 
summary( correct <- lmer(yij ~ xij + zj + (1 + xij |macro), data=R.intslope, REML=FALSE) )

# Too complex
summary( complex <- lmer(yij ~ xij + zj + xij*zj + (1 + xij  | macro), data=R.intslope, REML=FALSE) )

screenreg(list(simple,correct,complex))