---
title: 'Epds/Psyc/Stat 587: Inference for Random Effects'
author: "Carolyn J. Anderson"
date: '2023-02-06'
output:
  html_document:
    toc: yes
---
# Overview
This document mostly shows plots of random effects and how extract them.

# Setup

## Libraries
```{r libraries}
library(lme4)
library(lmerTest)
library(HLMdiag)
```

## Read in data
We will use data previously created so that we can skip creating 
various variables.

```{r indata}
setwd("D:\\Dropbox\\edps587\\lectures\\7 inference2")
#setwd("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\7 inference2")

hsb <- read.table("hsball.txt",header=TRUE)

# Number of schools will come in handy
N <- length(unique(hsb$id))
```

# Model we'll Work with

```{r}
summary( model.1 <- lmer (mathach ~ cSES + (1 |id),data=hsb,REML=FALSE) )
```

# Graph of EBLUPS by school +/- 1 sd
## Step 1: get setimates of random effects
Extract random effects, which yields a "mer"  object
```{r step1}
ranU <- ranef(model.1)
class(ranU)
```

## Step 2: change from mer to data.frame
Change class of "ranU" from mer to data.frame
```{r mer_to_data.frame}
df.1 <- as.data.frame(ranU)
head(df.1)
```

"df" is short for "data frame"

## Step 3: sort runU
Sort by ranU, which is "condval" (conditional value)
```{r sort_ranu}
df <- df.1[order(df.1$condval),]
```

## Step 4: integers to replace  school id
Add integers for school id's, which is needed to plot them in Step 5.
```{r int}
df$xgrp <- seq(1:N)
```

## Setp 5: Draw graph
```{r graph_ranu1}
plot(df$xgrp, df$condval, type="p",                      # plot points
     main="Model 1: Random intercepts +/- 1 std error",  # title
	   xlab="Schools (sort by Uo)",                       # label for horizontal axis
	   ylab="EBLUP of Uo",                                # label for vertical axis
	   col="blue")

arrows(df$xgrp, df$condval - df$condsd,       # origin
       df$xgrp, df$condval + df$condsd,       # terminus
	     length=0.00)
```

To draw in the sd errors, I drew arrows but set the arrowhead to size 0.  This is a hack solution, but it does the trick.

## Histogram of U_oj

```{r hist_ranU1}
hist(df.1$condval, 
     main="Estimated U_oj from Simple Model",
     col="gray",
     xlab="Emprical Bayes Estimates of U_oj")
```

# QQ plot of Us with 95% error bars                            
            

Simple QQplot for assessing normality of the $U_{0j}$s,
```{r qqplot}                          
qqdata <- qqnorm(df$condval)
```

A QQplot but with error bars.  This uses the same hack to get error bars in the plot.
```{r qqplot_with_error_bars}
plot(qqdata$x,qqdata$y,
     type="p", pch=19, col="red",
     main="Normal QQ plot of Uoj with 95% Confidence Bands",
     xlab="Theoretical Value",
     ylab="Estimated Uoj (EBLUPS)",
     ylim=c(-8,8))
     
arrows(qqdata$x, qqdata$y - 1.96*df$condsd, 
       qqdata$x, qqdata$y + 1.96*df$condsd, 
	     length=0.00, col="gray")

qqline(df$condval, col="steelblue")
```

# School Specific Regressions                                    

These fitted regressions for each school,

## Step 1:  Merge random effects and orginal data set
```{r}
df$id <- df$grp
hsb <- merge(hsb,df,by.x="id")
```

## Step 2:  Compute population average and school specific (conditional) values

### Population average values
```{r population_pop}
# This fixed effects
gamma <- as.data.frame(fixef(model.1))

hsb$g00 <- gamma[1,1]
hsb$g10 <- gamma[2,1]

# Population average model
hsb$ypop.ave     <- hsb$g00 + hsb$g10 * hsb$cSES
```

### School specific regressions
```{r school_specific}
# School specific (conditionally) regressions
hsb$yschool.spec <- hsb$g00 + hsb$g10 * hsb$cSES + hsb$condval

hsb$yschool.0 <- hsb$g00 +  hsb$condval
```

## step 3:  Plot population average

```{r Pop_average}
plot(hsb$cSES, hsb$ypop.ave, type='l',
     ylim=c(0,25),
     xlab="Centered SES",
  	 ylab="Fitted Math Scores",
	   main="Population Average model")
```

## step 4  
Plot separate regression line for each school.  These are the conditional regressions, conditional on the $\hat{U}_{0j}$s.
```{r school_plot_1}
# I create an empty plot to set up frame for graph
plot(hsb$cSES,hsb$ypop.ave,type='n',
     ylim=c(0,25),
     xlab="Centered SES",
	   ylab="Fitted Math Scores",
	   main="School Specific Regressions")
    hsb <- hsb[order(hsb$xgrp),]
    
# Now are in the subject specific regressions
for (i in (1:N)) {
   sub <- hsb[ which(hsb$xgrp==i),]
   lines(sub$cSES, sub$yschool.spec,type="l", col=i)
}
```    

Note that I used the index i to draw a different color for each regression


## step 4b:  Plot conditional regressions for 10 schools
Looking at all 160 regression lines is hard to see what's going on.  So we can take a subset of them to plot.
```{r just_10}
# Create empty frame 
plot(hsb$cSES,hsb$yschool.spec,type='n',
     ylim=c(0,25),
     xlab="Centered SES",
	   ylab="Fitted Math Scores",
	   main="School Specific Regressions, just 10 of them")

# Add in the 10 subgroups
hsb <- hsb[order(hsb$xgrp),]
for (i in (50:60)) {
   sub <- hsb[ which(hsb$xgrp==i),]
   lines(sub$cSES, sub$yschool.spec, type="l", col=i)
}
```

It would be better to take a random sample of schools, but in this case the schools are not ordered by any of the variables in the data set.


# QQPlot of population average values

```{r}
qqnorm(hsb$ypop.ave, pch=1, main="QQplot of Population Averages")

qqline(hsb$ypop.ave, col = "steelblue", lwd = 2)
```

# Histogram from HLMdiag

Get the residuals from HLMresid---this is not working even after changing depricated comdend. 
```{r hlmdiag, eval=FALSE}
res1 <- hlm_resid(model.1, level=1, type="EB", standardize=TRUE)
head(res1)
```

Now draw histogram and overlay a normal distribution.
```{r hist_normal1, eval=FALSE}
h<- hist(res1,breaks=15,density=20,main="Standardize Residuals")  # draw historgram
xfit <- seq(-40, 40, length=50)                  # sets range & number of quantiles
yfit <- dnorm(xfit, mean=0, sd=sd(res1))         # should be normal
yfit <- yfit*diff(h$mids[1:2])*length(res1)      # use mid-points 
lines(xfit, yfit, col='blue', lwd=2)             # draws normal
```

Alternative, use what we already found,
```{r hist_normal2}

## ----hist_normal2----------------------------------------------------
h<- hist((hsb$math-hsb$yschool.spec),
    breaks=15, col="lightblue1",
	main="Raw (conditional) Residuals: Rij")
```

Note that the values on x-axis are conditional residuals, conditional on $U_{0j}$s.

# Random Intercept and Slope Models

## Simple one
```{r RI_and_RS_model}
model.2 <- lmer (mathach ~ cSES + (1 +cSES |id),data=hsb,REML=FALSE)
summary(model.2)
```

## More Complex One

For the model that I wanted to fit, I had trouble with singularity.  To deal with this I multipled cSES by 1/10 to move the variance random effect for cSES away from 0. This worked.

To interpret effects in scale of cSES divide by 10

*         tau_cses/10 or tau_cses^2/100
*         gamma_cses/10

Clumsy but it works
```{r the_multiplication} 
hsb$mover <- rep(0.10, nrow(hsb))    
x <- as.matrix(hsb$cSES,  ncol=1)
y <- as.matrix(hsb$mover, ncol=1)
hsb$xyses <- x*y
```

Using re-scaled SES, fit model 3
```{r rescaled_ses_complex_model}
model.3 <- lmer(mathach ~ 1 + xyses + female + meanses + (1 + xyses + female | id),
             control = lmerControl(optimizer ="Nelder_Mead"), 
             data=hsb, REML=FALSE)
summary(model.3)
```

Extract random effects for this model
```{r extract_model3}
U <- as.data.frame(ranef(model.3))

Uoj <- U[which(U$term=="(Intercept)"),]
U1j <- U[which(U$term=="xyses"),]
U2j <- U[which(U$term=="female"),]
```

## QQplot U's from Complex Model

The reference lines were drawn in using a bit of trail and error

```{r qqplot_complex_model}
par(mfrow=c(2,2))          # a 2 x 2 matrix of plots

qqnorm(Uoj$condval,main="QQplot: Random Intercepts")
qqline(Uoj$condval,col="red")

qqnorm(U1j$condval,main="QQplot: Random cSES Effects")
qqline(U1j$condval,col="red")

qqnorm(U2j$condval,main="QQplot: Random female Effects")
qqline(U2j$condval,col="red")


x<- seq(0:100)
y<- seq(0:100)
plot(x,y, type="n", axes=FALSE, ylab="", xlab="")
text(55,100,"mathach ~ 1", cex=.9)
text(55,80," +cSES+female", cex=.9)
text(55,60,"+ meanses + (1",  cex=.9)
text(55,40,"+cSES+female|id)", cex=.9)
text(60,20,"Deviance=46494.5", cex=.9)
text(60,3,"AIC=46516.5", cex=.9)
#text(30,50,"tau_00 = 3.18")
#text(30.40,"Variance for cses...")
#text(30,30,"tau_11 = 1.49")
#text(30,20,
## ----qqplot_complex_model--------------------------------------------
par(mfrow=c(2,2))          # a 2 x 2 matrix of plots

qqnorm(Uoj$condval,main="Random Intercepts")
qqline(Uoj$condval,col="red")

qqnorm(U1j$condval,main="Random cSES Effects")
qqline(U1j$condval,col="red")

qqnorm(U2j$condval,main="Random female Effects")
qqline(U2j$condval,col="red")


x<- seq(0:100)
y<- seq(0:100)
plot(x,y, type="n", axes=FALSE, ylab="", xlab="")
text(55,100,"mathach ~ 1", cex=.95)
text(55,80," +cSES+female", cex=.95)
text(55,60,"+ meanses + (1",  cex=.95)
text(55,40,"+cSES+female|id)", cex=.95)
text(60,20,"Deviance=46494.5", cex=.95)
text(60,3,"AIC=46516.5", cex=.95)
#text(30,50,"tau_00 = 3.18")
#text(30.40,"Variance for cses...")
#text(30,30,"tau_11 = 1.49")
#text(30,20,"tau_22 = 0.91")
#text(30,10,"sigma^2=36.36")
```

I could not fit all text that I would have liked, but I did fit enough so that I would remember what model I'm looking at

## Bi-variate Scatter Plot

```{r}
# lattice plot of all random effects 
plot(ranef(model.3), main="Estimated Random Effects",col="blue")
```

# Shrinkage

The variance of $\hat{U}_{kj}$s are less than $\tau^{}$s.   The $\hat{U}_{kj}$ are shrunken toward the over all mean; that is, 
$$\hat{x}_{k+} \le |\hat{U}_{jk}| \le |\hat{x}_{kj}|$$.

The implication is the $$\mbox{var}(\hat{U}_{kj})  < \hat{\tau}_{k}^2$$

## As an illustration of shrinkage
Compare this to $\hat{\tau}_1^2= 8.612$ 
```{r shrinkage}
# model.1
var(df.1$condval)
```
Compare these to $\hat{\tau}_{k}^2= 3.1846$, $59.719$, and 0.9107$
```{r shrinkage_model_3}
#   model.3
var(Uoj$condval)
var(U1j$condval)
var(U2j$condval)
```

## Plot showing the shrinkage

This requires computing the grand mean for math, the school means for math, and the school specific math scores. The latter are conditional on the random effects. 

For this graphic, I will use model.1 and set cSES=0 (i.e., use yschool.0 here)

### Grand mean
```{r grand_mean}
hsb$grand <- rep(mean(hsb$mathach), nrow(hsb))
```

### Get mean math for each group
```{r compute_mean_math_for_each_school}
grp.mean <- aggregate(hsb$mathach, by=list(hsb$id), FUN=mean)
names(grp.mean) <- c("id","mean.math")
hsb <- merge(hsb, grp.mean, by="id")
```

### For a nice graph
The data need to be ordered; this will do nicely
```{r order_for_nice_graph}
hsb <- hsb[order(hsb$yschool.0),]
```

We also need an id for this ordering,
```{r new_id}
nj <- as.data.frame(table(hsb$id))
names(nj) <- c("id","count")
nj$integer <- 1:160
j <- 1
for (i in 1:160) {
  for (k in 1:nj$count[i]) {
     hsb$id.int[j] <- i
	 j <- j + 1
  }  
}	
```
## The graph
```{r shrinkage_graph}
par(mfrow=c(1,1))

plot(hsb$id.int, hsb$grand, type="l", col="steelblue",
     main="Shrinkgage: Model 1 when cses=0",
	 xlab="School ID",
	 ylab="Math",
	 cex=1.4,
	 lwd=2)

lines(hsb$id.int, hsb$yschool.0, col="red", type='l')
lines(hsb$id.int, hsb$mean.math, col="green",type="l")
legend(1, 18, legend=c("Grand Mean", "Conditional", "School Mean"),
       col=c("steelblue","red", "green"), lty=1:2, cex=1.1)
```
