---
title: "NELS 23 Schools and HLM"
author: "Carolyn J. Anderson"
date: "10/21/2022"
output: 
  html_document:
    toc: yes
---
Modifications:

*   For model 1, show how to add options to summary using add.summary

*   Re-did model 4 but center to show impact on effsize and time.


These examples use the NELS data from 23 (non-random) sample from the full 1003 schools in the data set (See Kreft & de Leeuw text).  lmer from the lme4 package will be used to fit model used REML and runjags will be used to fit Bayesian models.  I'm an only using MLE for comparison, which some of you are familar with. 

Some of these take a "long" time to run relative to the GIBBS we've done so far.

The variables in the data set are 

 Variable | Description           | Type    | Level
 ---------|-----------------------|---------|---------
SCHOOL    | SCHOOL ID             | integer | school
STUDENT   | STUDENT ID            | integer | student
SEX       | STUDENT SEX           | factor  | student 
RACE      | STUDENT RACE          | factor  | student
HOMEW     | TIME ON MATH HOMEWORK | numeric | student
          |                       | factor  | student
SCHTYPE   | SCHOOL TYPE           | factor  | student
SES       | SOCIOECONOMIC STATUS  | numeric | student                    
PARED     | PARENTAL EDUCATION    | factor  | student
          |                       | numeric | student
MATH      | MATH SCORE            | numeric | student
WHITE     | WHITE RACE BINARY     | factor  | student  
CLASSSTR  | CLASS STRUCTURE       | factor  | school
SCHSIZE   | SCHOOL SIZE           | numeric | school
URBAN     | URBANICITY            | factor  | school
GEO       | GEOGRAPHIC REGION     | factor  | school
MINORITY  | PERCENT MINORITY      | factor  | school
RATIO     | STUDENT-TEACHER RATIO | numeric | school
SCHTYPE   | School type           | factor  | school
private    | private SCHOOL BINARY  | factor  | school
------------------------------------------------------



# Set up

These are libraries that we'll use (at least most of them)
```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(lme4)     # MLE and REML of HLM
library(lmerTest) # test fixied effects  
library(coda)
library(rjags)
library(runjags)
library(jagsUI)   
library(mvtnorm)  # For correlated random effects
```

Set working directory and read in the data

```{r}
setwd("D:/Dropbox/edps 590BAY/Lectures/8 Multilevel models")

nels <- read.table("school23_data.txt",header=TRUE)
head(nels)
tail(nels)


nels$white[nels$race==4] <- 1
nels$white[nels$race<=3] <- 0
nels$white[nels$race>=5] <- 0

#
# ReCode schtype to private
#

nels$private <- NA
nels$private[nels$schtype==1] <- 0
nels$private[nels$schtype>=2] <- 1
```
## Various statistics

There are some that we need to set 

```{r}
# sort by schools:
nels <- nels[order(nels$school,nels$student) , ]

# Get information on number schools & create consecutive integer school.id
school <- unique(nels$school)                       #  
N <- length(school)                                 # need N and school.id
school.int <- as.data.frame(cbind(school,seq(1:N))) #  for dataList & model  
names(school.int) <- c("school", "school.id")       #
nels <- merge(nels,school.int,by="school")          #

# total number of students
n <- length(nels$math)                              # need n for dataList & model
```

```{r eval=FALSE}

## Double-check set-up
head(nels,n=50)
```

```{r}
# You don't need this but might want it
nj <- matrix(999,nrow=N,ncol=1)
for (j in 1:N) {
  nj[j] <- nrow(subset(nels,nels$school.id==j))
}  
```

## Descriptive Statistics

Should compute and examine to learn more about the variables.

## Graphing Data

We now have everything need to create a graph showing the variability between schools.

```{r}
# Graph to illustrate need for random intercept and slope
plot(nels$homew,nels$math,
	type="n", 
	col="blue",
	lwd=2,
	main="NELS: Linear Regression by School \n(for where there is data)",
	xlab="Time Spent Doing Homework",
	ylab="Math Scores"
	)
for (j in 1:N) {
  sub <- subset(nels,nels$school.id==j)
  fitted <- fitted(lm(math~homew,sub))
  lines(sub$homew,fitted,col=j)
}  
```

It looks like different intercepts and difference slopes for time spent doing homework.

# Model 1: Simple Random intercept, i.e., : math ~ 1 + homew + (1 |school.id)  

Now for our first model:

```{r}
###################################################################
# Random intercept     math~ 1 + homew + (1 |school.id)           #
###################################################################
#
# MLE
#
ri1.lmer <- lmer(math~ 1 + homew + (1 |school.id), data=nels, REML=TRUE)
```

We'll look at ri1.lmer later after we fit the model using Bayesian methods

```{r}
#
# Bayesian
#
N = 23

n = length(nels$math)

dataList <- list(
			y = nels$math,
			hmwk = nels$homew,
			school.id = nels$school.id,
			N    = N,
			n    = n,
			sdY  = sd(nels$math)
		 )
		 
ri.mod1 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)		 
		   mu[i] <-  g0 + U0j[school.id[i]] + g1*hmwk[i]  
		}		
	   for (j in 1:N) {
	     U0j[j] ~ dnorm(0,ptau)
	    }		
	    g0 ~ dnorm(0,1/(100*sdY^2))
	    g1 ~ dnorm(0,1/(100*sdY^2))

      tau ~ dunif(0.0001,200)
      ptau <- 1/tau^2
		
		  sigma ~ dunif(0.0001,2000)
      precision <- 1/sigma^2	
	}"

writeLines(ri.mod1, con="ri.mod1.txt")

start1 = list("g0"=mean(nels$math), "g1"=rnorm(1,1,3), "sigma"=sd(nels$math), "tau"=.5,   
			  .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,-2,3), "sigma"=runif(1,.001,10), "tau"=runif(1,0.0001,10),     
			  .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,0,3), "sigma"=runif(1,.001,10),  "tau"=runif(1,0.0001,10),     
			  .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=rnorm(1,-3,10), "g1"=rnorm(1,5,3), "sigma"=runif(1,.001,10),  "tau"=runif(1,0.0001,10),     
			  .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)
```

When writing the model, it is easy to make mistakes.  Rather than running code that could take a long time, a quick check of your code can save time.  The following is a fast check.

```{r eval=FALSE}
##############################################################
# Fast way to check your model:                              # 
#    rjags and then run for sampling jagsUI with parallel.   #
##############################################################
 
###################################################################################
```
Since I have already run this model and know the code is OK, I did not run the above.  So, let's use runjags (I'm lazy and runjags will compute statistics I want automatically, as well as able to do parallel).

```{r}
start.time <- Sys.time()
ri.mod1.runjags <- run.jags(model=ri.mod1,      
          method="parallel",  
          monitor=c("g0","g1","sigma", "tau"),
          data=dataList,
	        sample=10000,				  
		      n.chains=4,
		      inits=start)
end.time <- Sys.time()
(Ellapsed_time <- end.time - start.time)

plot(ri.mod1.runjags)
add.summary(ri.mod1.runjags,  confidence=c(.025, .975))
summary(ri1.lmer)		
```

Some graphics

```{r}
#-------------
# Graph to illustrate need for random intercept and slope
par(mfrow=c(2,2))
plot(nels$homew,nels$math,	type="n", lwd=2,
	main="NELS: Linear Regression by School \n(for where there is data)",
	xlab="Time Spent Doing Homework",	ylab="Math Scores"	)
for (j in 1:N) {
  sub <- subset(nels,nels$school.id==j)
  fitted <- fitted(lm(math~homew,sub))
  lines(sub$homew,fitted,col=j)
}  

# Compute conditional means (regressions)
S<-1
math.post <- matrix(999,nrow=S,ncol=n)
homework <- matrix(nels$homew,nrow=n,ncol=1)
for (s in 1:S) {
   pointer <- 1 
   for (j in 1:N) {
     U0 <-rnorm(1, 46.32, sd=4.9156)
     for (i in 1:nj[j]) {
       math.post[s,pointer] = U0 + 2.4006*homework[pointer] 
       pointer <- pointer + 1
   } 
  }
} 
nels$math.ri1 <- matrix(math.post,nrow=n,ncol=1)
plot(nels$homew,nels$math,	type="n", lwd=2,
	main="Fitted from Random Intercept \nmath~ 1 + homew + (1 |school.id)",
	xlab="Time Spent Doing Homework",	ylab="Math Scores"	)
for (j in 1:N) {
  sub <- subset(nels,nels$school.id==j)
  lines(sub$homew,sub$math.ri1,col=j)
}  
```

We knew that we needed more than just a random intercept and the graph also shows this...data (left) versus conditional fitted values (right).

# Model 2:  Add random Slope

So we add a random slope, but we'll start with uncorrelated random effects.

## Model 0:  Uncorrelated random effects
```{r}
###################################################################
# Random intercept & slope    math~ 1 + homew + (1 |school.id)    #
#  Un-correlated random effects                                   #
###################################################################
ris1.lmer <- lmer(math~ 1 + homew + (1 |school.id) + (0 + homew | school.id), data=nels, REML=TRUE)
summary(ris1.lmer)

dataList <- list(
			math = nels$math,
			hmwk = nels$homew,
			school.id = nels$school.id,
			N    = N,
			n    = n,
			sdY  = sd(nels$math)
		 )
		 
ris.mod1 <- "model { 
       for (i in 1:n) {
	     math[i] ~ dnorm(mu[i],precision)
		 mu[i] <- g0 + g1*hmwk[i] + b0j[school.id[i]] + b1j[school.id[i]]*hmwk[i]
		}
		
	   for (j in 1:N) {
	     b0j[j] ~ dnorm(0,ptau0)
	     b1j[j] ~ dnorm(0,ptau1)
	    }
		
	    g0 ~ dnorm(0,1/(100*sdY^2))
	    g1 ~ dnorm(0,1/(100*sdY^2))

        tau0 ~ dunif(0.0001,200)
        ptau0 <- 1/tau0^2
		
        tau1 ~ dunif(0.0001,200)
        ptau1 <- 1/tau1^2
		
		sigma ~ dunif(0.0001,2000)
        precision <- 1/sigma^2	
		}"

writeLines(ris.mod1, con="ris.mod1.txt")

start1 = list("g0"=mean(nels$math), "g1"=rnorm(1,1,3), 
              "sigma"=sd(nels$math),  "tau0"=.5, "tau1"=.5,  
			  .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,-2,3), 
              "sigma"=runif(1,.001,10), "tau0"=runif(1,0.0001,10), "tau1"=runif(1,0.0001,10),     
			  .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,0,3),
              "sigma"=runif(1,.001,10),  "tau0"=runif(1,0.0001,10), "tau1"=runif(1,0.0001,10),         
			  .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=rnorm(1,-3,10), "g1"=rnorm(1,5,3), 
              "sigma"=runif(1,.001,10),  "tau0"=runif(1,0.0001,10), "tau1"=runif(1,0.0001,10),     
			  .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)
```

```{r eval=FALSE}
###########################################################################
# Fast way to check your model:  rjags and then run for sampling runjags  #
#   with parallel.                                                        #
###########################################################################
ris.mod1 <- jags.model(file="ris.mod1.txt", # compiles and initializes model
                        data=dataList,
                        inits=start1,
                        n.chains=4,
                        n.adapt=500) 
 
###################################################################################
```

We'll go ahead and run this,

```{r}
start_time <- Sys.time()
ris.mod1 <- run.jags(model=ris.mod1,      
          method="parallel",  
          monitor=c("g0","g1","sigma", "tau0","tau1"),
          data=dataList,
          sample=10000,				  
		      n.chains=4,
		      inits=start,
	        thin=10)
end_time <- Sys.time()
(ellapsed.time <- end_time - start_time)


summary(ris.mod1)
```

And for comparison
```{r}
summary(ris1.lmer)
```
And some graphs from Bayesian analysis
```{r, eval=FALSE}
plot(ris.mod1.runjags)
```

## Mone Carlo to obtain estmates of the random effects

Note coefficients used should match what you just got...
```{r} 
#### Better to estimate U within jags, but to get idea across ######
# Compute conditional regressions
S<-1
math.post <- matrix(999,nrow=S,ncol=n)
homework <- matrix(nels$homew,nrow=n,ncol=1)
for (s in 1:S) {
   pointer <- 1 
   for (j in 1:N) {
     U0 <-rnorm(1, 46.443, sd=7.5859)
	   U1 <-rnorm(1, 1.9677, sd=3.9674)
     for (i in 1:nj[j]) {  # for estimates include + rnorm(1,g0,sigma)
       math.post[s,pointer] = U0 + U1*homework[pointer]   
       pointer <- pointer + 1
   } 
  }
} 
```

Now we can plot stuff

```{r}
nels$math.ris12 <- matrix(math.post,nrow=n,ncol=1)
plot(nels$homew,nels$math,	type="n", 	lwd=2,
	main="Fitted from Random Intercept & Slope \nmath~1+homew+(1|school.id)+(0+homew|school.id)",
	xlab="Time Spent Doing Homework",	ylab="Math Scores",
	ylim=c(30,70)	)
for (j in 1:N) {
  sub <- subset(nels,nels$school.id==j)
  lines(sub$homew,sub$math.ris1,col=j)
}  
```

This look a lot like the data, but we can do better by allowing random intercept and random slope to be correlated.  

# Model 2: correlated random effects
```{r}

#############################################################
# Random intercept & slope                           Wishart#
#   math~ 1 + homew + (1 + homew|school.id)                 #
#   with correlated random effects:  use wishart            #
#############################################################
ris2.lmer <- lmer(math~ 1 + homew + (1 + homew|school.id), data=nels, REML=TRUE)
summary(ris2.lmer)

dataList <- list(
			y = nels$math,
			hmwk = nels$homew,
			school.id=nels$school.id,
			n    = n,
			N    = N,
			sdY  = sd(nels$math)
		 )
		 
re.mod2 <- "model { 
    # Likelihood: the data model
        for (i in 1:n) {
	       y[i] ~ dnorm(meanY[i],precision)
		     meanY[i] <-  betaj[school.id[i],1] + betaj[school.id[i],2]*hmwk[i]
		   }
	# Random Effects -- note multivariate normal density is used
	   for (j in 1:N) {
	     betaj[j,1:2] ~ dmnorm(mu[1:2],Omega[1:2,1:2])
	    }
		
	# Priors
     precision ~ dgamma(0.01,0.01)
	   sigma <- 1/sqrt(precision)	
	   
	   mu[1] ~ dnorm(0,1/(100*sdY^2))# or mu[i]<- 0 and add g0 & g1 to meanY 
	   mu[2] ~ dnorm(0,1/(100*sdY^2))# and add g0<- dnorm(0,1,(100*sdY^2))
 	                                # & for g1.   This will make it easier
		# to add models for intercept & slopes.
		# See next model.
       Omega[1:2,1:2] ~ dwish(R[,],2.1)
       R[1,1] <- 1/2.1
       R[1,2] <- 0         # Un-correlated--there are alternative
       R[2,1] <- 0         # ways of setting this up by just adding
       R[2,2] <- 1/2.1     # uni-variate distribution for each tau
	     Tau <- inverse(Omega) # This turn into our taus.
		}"

writeLines(re.mod2, con="re.mod2.txt")

start1 = list("precision"=1, .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("precision"=1, .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("precision"=1, .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("precision"=1, .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

##########################################################################
# Fast way to check your model:  
# rjags and then run for sampling runjags with parallel.
##########################################################################
re.mod2x <- jags.model(file="re.mod2.txt",     # compiles and initializes model
                        data=dataList,
                        inits=start1,
                        n.chains=4,
                        n.adapt=500)            
#########################################################################
start_time <- Sys.time()
re.mod2.runjags <- run.jags(model=re.mod2,      
                  method="parallel",  
                  monitor=c("mu","sigma","Tau"),
                  data=dataList,
				  sample=20000,				  
				  n.chains=4,
				  inits=start)
end_time <- Sys.time()
(ellapsed.time <- end_time - start_time)


add.summary(re.mod2.runjags)
plot(re.mod2.runjags)
summary(ris2.lmer)
```		 

The frequentist estimated more 'failed to converge'.   It doesn't look that bad and we 
might be able to get it to converge if we change the optimization method.

Let's look at what fitted values look like conditional on random effects

```{r}

## Some Monte Carlo

############################################################ 
# Better to estimate U within jags, but to get idea across #
# Compute conditional regressions                          #
############################################################
S<- 1
math.post <- matrix(999,nrow=S,ncol=n)
homework <- matrix(nels$homew,nrow=n,ncol=1)
Tau <- matrix(c(63.01,-27.717,-27.717,17.804),nrow=2,ncol=2)
for (s in 1:S) {
   pointer <- 1 
   for (j in 1:N) {
    Uj <- rmvnorm(1,mean=c(46.296,2.0062),sigma=Tau) 
     for (i in 1:nj[j]) {
       math.post[s,pointer] = Uj[1,1] + Uj[1,2]*homework[pointer]   # for estimates include + rnorm(1,g0,sigma)
       pointer <- pointer + 1
   } 
  }
} 

nels$math.ris2 <- matrix(math.post,nrow=n,ncol=1)

plot(nels$homew,nels$math,
	type="n", 
	col="blue",
	lwd=2,
	main="Fitted from Random Intercept \nmath~ 1 + homew + (1 + homew|school.id)",
	xlab="Time Spent Doing Homework",
	ylab="Math Scores",
	ylim=c(30,70)
	)
for (j in 1:N) {
  sub <- subset(nels,nels$school.id==j)
  lines(sub$homew,sub$math.ris2,col=j)
}  		 
```

# Model 3: Cross-level Interaction

Now for model 3:
```{r}
#############################################################################
# Random intercept & slope with more predictors including a cross-level     #
#  interaction                                                              #
#    math~ 1 + homew + ses + private + homew*private +(1 + homew|school.id)   #
#                                                                           #
#############################################################################
# Turn into dummy codes
nels$private <- ifelse(nels$schtype==1,1,0)

ris3.lmer <- lmer(math~ 1 + homew + ses + private  + homew*private 
                  + (1 + homew|school.id), data=nels, REML=TRUE)
summary(ris3.lmer)

dataList <- list(
			math = nels$math,
			hmwk = nels$homew,
			ses  = nels$ses,
			private=nels$private,
			school.id=nels$school.id,
			n    = n,
			N    = N,
			sdY  = sd(nels$math)
		 ) 
	
re.mod3 <- "model { 
    # Likelihood: the data model
        for (i in 1:n) {
	       math[i] ~ dnorm(mmath[i],precision)
		   mmath[i] <- g0 +  Uj[school.id[i],1] + g3*private[i] 
		             + (g10 + Uj[school.id[i],2] + g11*private[i])*hmwk[i] + g2*ses[i] 
		   }  
	# Random Effects -- note multivariate normal density is used
	   for (j in 1:N) {
	    Uj[j,1:2] ~ dmnorm(mu[1:2],Omega[1:2,1:2])
 	    }
  	
	# Priors
     precision ~ dgamma(0.01,0.01)
	   sigma <- 1/sqrt(precision)
	   
  	 g0 ~ dnorm(0,1/(100*sdY^2))
   	 g2 ~ dnorm(0,1/(100*sdY^2))
	   g3 ~ dnorm(0,1/(100*sdY^2))
	   g10 ~ dnorm(0,1/(100*sdY^2))
	   g11 ~ dnorm(0,1/(100*sdY^2))
   
	   mu[1] <- 0     # This identifies the model 
     mu[2] <- 0     # 
	   
     Omega[1:2,1:2] ~ dwish(R[,],2.1)
       R[1,1] <- 1/2.1
       R[1,2] <- 0
       R[2,1] <- 0
       R[2,2] <- 1/2.1
	   Tau <- inverse(Omega)
		}"

writeLines(re.mod3, con="re.mod3.txt")

start1 = list("precision"=1,   "g0"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10), .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
start2 = list("precision"=0.5, "g0"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10), .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("precision"=0.05,"g0"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10),.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("precision"=.002,"g0"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10),.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
```

If you want to check whether model is OK, then run
```{r eval=FALSE}
##################################################################################
# Fast way to check your model:   rjags and then run for sampling runjags with parallel.
 re.mod3 <- jags.model(file="re.mod3.txt",     # compiles and initializes model
                        data=dataList,
                        inits=start1,
                        n.chains=4,
                        n.adapt=500)            
###################################################################################
```

If model is OK, then run jags
```{r}
start_time <- Sys.time()
re.mod3.runjags <- run.jags(model=re.mod3,      
           method="parallel",  
           monitor=c("g0", "g10", "g11","g2","g3","sigma","Tau"),
           data=dataList,
				   sample=10000,				  
			    n.chains=4,
				  inits=start,
				  thin=10)
end_time <- Sys.time()
(ellapsed_time <- end_time - start_time)
#
# run these in class but not for Rmd 
plot(re.mod3.runjags)    
add.summary(re.mod3.runjags)

summary(ris3.lmer)
```

# Model evaluation: Model 4


I will just put in large chunks of code and we'll discuss in class what they do.

```{r}
dataList <- list( 	math = nels$math,
			hmwk = nels$homew,
			ses  = nels$ses,
			private=nels$private,
			school.id=nels$school.id,
			n    = n,
			N    = N,
			sdY  = sd(nels$math)
		 ) 
	
re.mod4 <- "model { 
    # Likelihood: the data model
        for (i in 1:n) {
	       math[i] ~ dnorm(mmath[i],precision)
		     mmath[i] <- g00 +  Uj[school.id[i],1] + g3*private[i] 
		       + (g10 + Uj[school.id[i],2] + g11*private[i])*hmwk[i] + g2*ses[i] 

		   intercept.blup[i] <-  Uj[school.id[i],1]    # intercept random effect(posterior)
		   slope.blup[i]     <-  Uj[school.id[i],2]    # slope random effect
		   marg.post[i]      ~ dnorm(mmath[i],precision) 
		}  
	# Random Effects -- note multivariate normal density is used
	   for (j in 1:N) {
	    Uj[j,1:2] ~ dmnorm(mu[1:2],Omega[1:2,1:2])
 	    }
  	
	# Priors
     precision ~ dgamma(0.01,0.01)
	   sigma <- 1/sqrt(precision)
	   
  	 g00 ~ dnorm(0,1/(100*sdY^2))
   	 g2 ~ dnorm(0,1/(100*sdY^2))
	   g3 ~ dnorm(0,1/(100*sdY^2))
	   g10 ~ dnorm(0,1/(100*sdY^2))
	   g11 ~ dnorm(0,1/(100*sdY^2))
   
	   mu[1] <- 0     # This identifies the model 
     mu[2] <- 0     # 
	   
     Omega[1:2,1:2] ~ dwish(R[,],2.1)
       R[1,1] <- 1/2.1
       R[1,2] <- 0
       R[2,1] <- 0
       R[2,2] <- 1/2.1
	   Tau <- inverse(Omega)
		}"

writeLines(re.mod4, con="re.mod4.txt")

start1 = list("precision"=1,   "g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10), .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
start2 = list("precision"=0.5, "g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10), .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("precision"=0.05,"g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10),.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("precision"=.002,"g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10),.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
```

Now get samples but monitor main effects (otherwise I have 1567 variables x 15000 samples).

```{r}
start_time <- Sys.time()
re.mod4.runjags <- run.jags(model=re.mod4,      
                  method="parallel",  
                  monitor=c("g00", "g10", "g11","g2","g3","sigma","Tau"),
                  data=dataList,
				  sample=15000,				  
				  n.chains=4,
				  inits=start,
				  thin=10)
end_time <- Sys.time()
(ellapsed.time <- end_time - start_time)

plot(re.mod4.runjags)
gelman.plot(re.mod4.runjags)
geweke.diag(re.mod4.runjags, frac1=.10, frac2=.50)
add.summary(re.mod4.runjags)  
```

Now monitor all
```{r}
start_time <- Sys.time()
re.mod4.runjags <- run.jags(model=re.mod4,      
                  method="parallel",  
                  monitor=c("g00", "g10", "g11", "g2", "g3", "sigma", "Tau", "marg.post", "intercept.blup", "slope.blup"),
                  data=dataList,
				  sample=15000,				  
				  n.chains=4,
				  inits=start,
				  thin=10)
end_time <- Sys.time()
(ellpased_time <- end_time - start_time)

add.summary(re.mod4.runjags)  
```
Now we can look at posterior starting with some data manipulation

```{r}
ffit <- mcmc(re.mod4.runjags)
names(ffit)                    # lets you know what is available

# posterior samples are in ffit$mcmc ---> there are separate columns for each chain, i.e., 15000 x 4
g00 <- as.array(ffit$mcmc[,1])
g10 <- as.array(ffit$mcmc[,2])
g11 <- as.array(ffit$mcmc[,3])
g2  <- as.array(ffit$mcmc[,4])
g3  <- as.array(ffit$mcmc[,5])
sigma <- as.array(ffit$mcmc[,6])
tau11 <- as.array(ffit$mcmc[,7])
tau12 <- as.array(ffit$mcmc[,8])
tau22 <- as.array(ffit$mcmc[,10])
```

Some more data manipulations

```{r}
margin.post <- as.array(ffit$mcmc[,11:529])       # separate matrices for each chain, i.e. 10000 x 519 x 4
intercept.blup <- as.array(ffit$mcmc[,530:1048])   
slope.blup <- as.array(ffit$mcmc[,1049:1567])  
                                                  # just want fit and fit.new
# to get 40,000 x 519: 
marginal.means <- rbind(margin.post[,1:519,1],margin.post[,1:519,2],margin.post[,1:519,3],margin.post[,1:519,4])
int.blup       <- rbind(intercept.blup[,1:519,1],intercept.blup[,1:519,2],intercept.blup[,1:519,3],intercept.blup[,1:519,4])
slp.blup       <- rbind(slope.blup[,1:519,1],slope.blup[,1:519,2],slope.blup[,1:519,3],slope.blup[,1:519,4])


# Sum over rows to get fitted values (marginal) for each student
nels$student.margin <- apply(marginal.means,2,"mean")				  
```

Using computed values above, plot data and fitted marginal and conditional distributions

```{r}
# Data & fitted margin
par(mfrow=c(2,2))
hist(nels$math,breaks=20,main="Data: Math margin")
hist(nels$student.margin,breaks=20,main="Marginal Distribution")

   											  
# to get the conditional margin: 
cond.mean <- marginal.means + int.blup + slp.blup
nels$conditional.means <- apply(cond.mean,2,"mean")												  
												  
hist(nels$conditional.mean,breaks=20,main="Conditional Means")												  
cor(nels$conditional.mean,nels$math)
```

Let's go the other way: collapse over students

```{r}
# Compute statistics over columns (students)
ymeans <- apply(cond.mean,1,"mean")
ymed   <- apply(cond.mean,1,"median")
ymax  <- apply(cond.mean,1,"max")
ysd   <- apply(cond.mean,1,"sd")

# Bayesian p-values
pmean <- mean(mean(nels$math) < ymeans)
pmed  <- mean(median(nels$math) < ymed)
pmax  <- mean(max(nels$math) < ymax)
psd   <- mean(sd(nels$math) < ysd)
```

And some graphs for these

```{r}
# Graphs of descriptive statistics from posterior predictive distribution 
par(mfrow=c(2,2))
 
hist(ymed,breaks=10,main=paste("Conditional Medians", "\nBayesian P-value =", round(pmed, 2)))
lines(c(median(nels$math),median(nels$math)),c(0,10000),col="blue",lwd=2)

hist(ymax,breaks=10,main=paste("Conditional  Maximums", "\nBayesian P-value =", round(pmax, 2)))
lines(c(max(nels$math),max(nels$math)),c(0,10000),col="blue",lwd=2)

hist(ymeans,breaks=10,main=paste("Conditional  Means", "\nBayesian P-value =", round(pmean, 2)))
lines(c(mean(nels$math),mean(nels$math)),c(0,10000),col="blue",lwd=2)

hist(ysd,breaks=10,main=paste("Conditional  SDs", "\nBayesian P-value =", round(psd, 2)))
lines(c(sd(nels$math),sd(nels$math)),c(0,10000),col="blue",lwd=2)
```

And if you want to look at the estimated random effects

```{r}
plot(nels$math,nels$conditional.means,type="p")
plot(nels$math,nels$student.margin,type="p")

hist(int.blup)
hist(slp.blup)
```

# Re-do Model 4 but with Centering

```{r}
hmwk = nels$homew - mean(nels$homew)
ses  = nels$ses - mean(nels$ses)
private= nels$private - mean(nels$private)
N = 23
n = length(nels$math)

dataList <- list( 	math = nels$math,
			hmwk = nels$homew ,
			ses  = nels$ses ,
			private=nels$private,
			school.id=nels$school.id,
			n    = n,
			N    = N,
			sdY  = sd(nels$math)
		 ) 
	
re.mod4 <- "model { 
    # Likelihood: the data model
        for (i in 1:n) {
	       math[i] ~ dnorm(mmath[i],precision)
		     mmath[i] <- g00 +  Uj[school.id[i],1] + g3*private[i] 
		       + (g10 + Uj[school.id[i],2] + g11*private[i])*hmwk[i] + g2*ses[i] 

		   intercept.blup[i] <-  Uj[school.id[i],1]    # intercept random effect(posterior)
		   slope.blup[i]     <-  Uj[school.id[i],2]    # slope random effect
		   marg.post[i]      ~ dnorm(mmath[i],precision) 
		}  
	# Random Effects -- note multivariate normal density is used
	   for (j in 1:N) {
	    Uj[j,1:2] ~ dmnorm(mu[1:2],Omega[1:2,1:2])
 	    }
  	
	# Priors
     precision ~ dgamma(0.01,0.01)
	   sigma <- 1/sqrt(precision)
	   
  	 g00 ~ dnorm(0,1/(100*sdY^2))
   	 g2 ~ dnorm(0,1/(100*sdY^2))
	   g3 ~ dnorm(0,1/(100*sdY^2))
	   g10 ~ dnorm(0,1/(100*sdY^2))
	   g11 ~ dnorm(0,1/(100*sdY^2))
   
	   mu[1] <- 0     # This identifies the model 
     mu[2] <- 0     # 
	   
     Omega[1:2,1:2] ~ dwish(R[,],2.1)
       R[1,1] <- 1/2.1
       R[1,2] <- 0
       R[2,1] <- 0
       R[2,2] <- 1/2.1
	   Tau <- inverse(Omega)
		}"

writeLines(re.mod4, con="re.mod4.txt")

start1 = list("precision"=1,   "g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10), .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
start2 = list("precision"=0.5, "g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10), .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("precision"=0.05,"g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10),.RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("precision"=.002,"g00"=rnorm(1,0,10),"g2"=rnorm(1,0,10), "g3"=rnorm(1,0,10),"g10"=rnorm(1,0,10), "g11"=rnorm(1,0,10),.RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
start <- list(start1,start2,start3,start4)
```

Now get samples but monitor main effects (otherwise I have 1567 variables x 15000 samples).

```{r}
start_time <- Sys.time()
re.mod4.runjags <- run.jags(model=re.mod4,      
                  method="parallel",  
                  monitor=c("g00", "g10", "g11","g2","g3","sigma","Tau"),
                  data=dataList,
				  sample=15000,				  
				  n.chains=4,
				  inits=start,
				  thin=10)
(end_time <- Sys.time())
(ellapsed_time <- end_time - start_time)

plot(re.mod4.runjags)
gelman.plot(re.mod4.runjags)
geweke.diag(re.mod4.runjags, frac1=.10, frac2=.50)
add.summary(re.mod4.runjags)  
```

Now monitor all
```{r}
start_time <- Sys.time()
re.mod4.runjags <- run.jags(model=re.mod4,      
                  method="parallel",  
                  monitor=c("g00", "g10", "g11", "g2", "g3", "sigma", "Tau", "marg.post", "intercept.blup", "slope.blup"),
                  data=dataList,
				  sample=15000,				  
				  n.chains=4,
				  inits=start,
				  thin=10)
end_time <- Sys.time() 
(ellpased_time <- end_time - start_time)
  
add.summary(re.mod4.runjags)  
```
Now we can look at posterior starting with some data manipulation