---
title: 'Edpsy/Psych/Stat 587:  Modeling Building'
author: "Carolyn J. Anderson"
date: "3/13/2023"
output: 
     html_document:
       toc: yes
---
We will start with exploratory data analysis, followed by model the data, and assessing model fit and assumptions.   Included in this example are location-scale models.

Many of the figures are in the lecture notes on model building.  The page numbers are those in the lecture notes (at least a close approximation of their location)

# Set up


```{r load_libraries, message=FALSE}
library(lme4)
library(lmerTest)
library(lattice)
library(ggplot2)
library(stringi)

library(HLMdiag)
library(texreg)
library(optimx)
```

We will use some of the functions described in lecture on inference of the marginal model (e.g., robust, contrast, hlmRsq, etc.)

```{r source_location_set}
source("D:\\Dropbox\\edps587\\lectures\\6 inference1\\All_functions.txt")
#source("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\6 inference1\\All_functions.txt")
```
# NELS Data

```{r read_in_data}
nels<-read.table("D:\\Dropbox\\edps587\\lectures\\8 modelbuilding\\nels23_data.txt",header=TRUE)
#nels<-read.table("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\8 modelbuilding\\nels23_data.txt",header=TRUE)

#look at first few rows of data set
head(nels)
```

Read in the data and create some variables we'll need.

```{r set_up_variables}
nels$public <- ifelse (nels$schtype==1, 1, 0)
nels$white  <- ifelse (nels$race==4, 1, 0)


schMses <- as.data.frame(aggregate(ses~school, data=nels, FUN="mean"))
names(schMses) <- c('school', 'schMses')
nels <- merge(nels,schMses, by=c('school'))
nels$schCses <- nels$ses - nels$schMses

nels$homewrk <- as.numeric(nels$homew)
nels$school  <- as.factor(nels$school)

head(nels)
```

That should be enough for now.

# Math as a function of Homework

## Linear regression model and overall 95% CI fit to the data
```{r overall_fit_single_level}
math.homewrk <- lm(math~homewrk, data=nels)
conflm1<-confint(math.homewrk)
plot(nels$homewrk,nels$math,type='p', main='Overall regression with 95% CI')
abline(math.homewrk)
abline(coef=conflm1[,1],lty=2)
abline(coef=conflm1[,2],lty=2)
```


##  Smooth (loess) curve fit to data
  
  Fit loess to data are plot
```{r fit_loess1}
lw1 <- loess(math ~ homewrk,data=nels)
# --- plot data
plot(nels$homewrk,nels$math,
	 main="NELS:Math Score Plotted vs Homework \n with smooth (loess) curve",
	 xlab="Time Spent Doing homework",
	 ylab="Math Scores",
	 cex=1.0,
	 col="black"
	 )	
# --- put in proper order
	j <- order(nels$homewrk)
# --- Add loess to figure	
lines(nels$homewrk[j],lw1$fitted[j],col="blue",lwd=3)
```


## Joint points for each school 
 This graph is around page 25.  We compute means mean of math for each level of time spent doing home for rach school.  These are plotted just connecting points (one curve for each school)

```{r means_by_school}
math.hmwk <- aggregate(math ~ homewrk + school, data=nels, FUN="mean")

ggplot(math.hmwk, aes(x=homewrk, y=math, col=school,type='l'),
  main='Regressions all in one plot')  + geom_line() +  theme(legend.position="none")
```

Or an alternative method

```{r alternative_method}
plot(nels$homewrk,nels$math,
	 main="NELS:Math Score Plotted vs Homework \n joining means for each school",
	 xlab="Time Spent Doing homework",
	 ylab="School Mean Math Scores",
	 cex=1.0,
	 col="black",
	 type="n"
	 )	
# --- put in proper order
	j <- order(math.hmwk$homewrk)
# --- Add to figure	
for (i in unique(math.hmwk$school)) {
  sub <- math.hmwk[which(math.hmwk$school==i),]
  lines(sub$homewrk,sub$math,col=i)
} 
```


# Lattice Plots 

## page 32 linear regression math~homewrk + race

```{r lattice_1}
xyplot(math ~ homew | school, data=nels, col.line='black', 
                      type=c('p'),
                      main='Varability in Math ~ homework relationship')
```

## page 19 - 20 
 Lattice plot with regressions overlaid

```{r lattice_2}
xyplot(math ~ homew | school, data=nels, col.line='black', 
                      type=c('p','r'),
                      main='Varability in Math ~ homework relationship')
```


## page 23-24   
Lattice plot with smooth regressions  (i.e., loess)

```{r lattice_3, warning=FALSE, message=FALSE}
xyplot(math ~ homew | school, data=nels, col.line='black', 
                      type=c('p','smooth'),
                      main='Varability in Math ~ homework relationship')
```

## page 25  
compute mean and plot by school just connecting points
```{r mean_math_join}
math.hmwk <- aggregate(math ~ homewrk + school, data=nels, FUN="mean")

ggplot(math.hmwk, aes(x=homewrk, y=math, col=school,type='l'),
  main='Regressions all in one plot')  + geom_line() +  theme(legend.position="none")
```
## page 26  All regressions in one figure

```{r overlay_regressions_hmwk}
ggplot(nels, aes(x=homewrk, y=math, col=school,type='l'),
       main='Regressions all in one plot')  + geom_smooth(method="lm", se=FALSE ) +  theme(legend.position="none")
```

## Alternate code to put all regressions in a single plot:
 Set up graphic frame:
```{r alternate_overlay}
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"
	)
# --- loop through the schools and add line
for (j in unique(nels$school)) {
  sub <- nels[ which(nels$school==j),]
  fitted <- fitted(lm(math~homew,sub))
  lines(sub$homew,fitted,col=j)
}   
```


# Math ~ homework + white
Need this for next set of graphs,

Draw mean math and plot by white
```{r get_means_by_white}
w <- as.factor(nels$white)
math.white <- aggregate(math ~ homewrk + w, data=nels, FUN="mean")
```

```{r graph_them}
ggplot(math.white, aes(x=homewrk, y=math, col=w, type='l'),
            main='Regressions all in one plot')  + geom_line() 
```

Alternate that looks nicer:
```{r page28_nicer}
plot(math.white$homewrk[1:7],math.white$math[1:7],type="l",
     main="Mean Math by Homework",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score",
	 col="blue")
lines(math.white$homewrk[8:14],math.white$math[8:14],type="l",col="red")
legend("topleft",legend=c("Non-White", "White"),
       col=c("blue","red"), lty=1, cex=1.2)   
```


## page 29 linear regression math~homewrk + white

```{r page29}
plot(math.white$homewrk,math.white$math,type="n",
     main="Regression of Math on Homework",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score",
	 ylim=c(40,70))
subw <- math.white[which(math.white$w==0),]
lm1 <- lm(math ~ homewrk,data=subw)
lines(subw$homewrk,fitted(lm1),col="red")
subw <- math.white[which(math.white$w==1),]
lm2 <- lm(math ~ homewrk,data=subw)
lines(subw$homewrk,fitted(lm2),col="blue")
legend("topleft",legend=c("Non-White", "White"),
       col=c("red","blue"), lty=1, cex=1.2)   
``` 

## Regression to whole uncollapsed data set

```{r yet_another_one}
plot(math.white$homewrk,math.white$math,type="n",
     main="Regression of Math on Homework",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score",
	 ylim=c(40,70))
subw <- nels[which(nels$white==1),]
abline(lm(math~homewrk,data=subw),col="blue")
subn <- nels[which(nels$white==0),]
abline(lm(math~homewrk,data=subn),col="red")
legend("topleft",legend=c("Non-White", "White"),
       col=c("red","blue"), lty=1, cex=1.2)   
```


#  Math mean ~ homework +  by race

Need to get means by race
```{r math_means_by_race}
r <- as.factor(nels$race)
math.race <- aggregate(math ~ homewrk + r, data=nels, FUN="mean")
```

## math ~ homework + race (page 32ish)

```{r plot_alternate_page32}
plot(math.race$homewrk[1:5],math.race$math[1:5],type="l",
     main="Mean Math by Homework",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score",
	 ylim=c(30,70),
	 xlim=c(0,7),
	 col="blue")
lines(math.race$homewrk[6:11],math.race$math[6:11],type="l", col="red")
lines(math.race$homewrk[12:18],math.race$math[12:18],type="l", col="darkgreen")
lines(math.race$homewrk[19:26],math.race$math[19:26],type="l", col="purple")
lines(math.race$homewrk[27:28],math.race$math[27:28],type="l", col="cyan")
legend("topleft",legend=c("Asian/PI", "Hispanic","Black","White","Native Amercian"),
       col=c("blue","red","darkgreen","purple","cyan"), lty=1, cex=1.2)   
```

## Some more

page 32 linear regression math~homework + race

regressions extend beyond range where there is data
```{r plot_32, warning=FALSE, message=FALSE}
par(mfrow=c(1,1))
plot(math.race$homewrk[1:5],math.race$math[1:5],type="n",
     main="Mean Math by Homework",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score",
	 ylim=c(30,70),
	 xlim=c(0,7))
abline(lm(math.race$math[1:5]~math.race$homewrk[1:5]),type="l",col="blue")
abline(lm(math.race$math[6:11]~math.race$homewrk[6:11]),type="l",col="red")
abline(lm(math.race$math[12:18]~math.race$homewrk[12:18]),type="l",col="darkgreen")
abline(lm(math.race$math[19:26]~math.race$homewrk[19:26]),type="l",col="purple")
abline(lm(math.race$math[27:28]~math.race$homewrk[27:28]),type="l",col="cyan")
legend("topleft",legend=c("Asian/PI", "Hispanic","Black","White","Native American"),
       col=c("blue","red","darkgreen","purple","cyan"), lty=1, cex=0.9)   
```

## page 34  race x homework  to show where there is not much data


Does not extend beyond the range of data.  Set up graphic frame and then add stuff
```{r page34}
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"
	)
# --- loop through the schools and add line
for (j in nels$race) {
  sub <- nels[ which(nels$race==j),]
  fitted <- fitted(lm(math~homew,sub))
  lines(sub$homew,fitted,col=j)
}
legend("topleft",legend=c("Asian/PI", "Hispanic","Black","White","Native Amercian"),
       col=c("blue","red","darkgreen","purple","cyan"), lty=1, cex=0.9)   
```
  
  
##  Cross-classification of homew x race

There's not much data for some races and homework levels
```{r cross-classification}
with(nels, table(homew, race))
```
 

# Math ~ homework + gender
  
  
##  page 36 mean math by gender
```{r mean_gender}
nels$gender <- as.factor(nels$sex)
math.gender <- aggregate(math ~ homewrk + gender, data=nels, FUN="mean")

nels$gender <- as.factor(nels$sex)
math.gender <- aggregate(math ~ homewrk + gender, data=nels, FUN="mean")
```

For many of these graphs, I also have code for ggplot; however, it doesn't seem to be 
working, for example...
```{r bad_ggplot, eval=FALSE}
nels$gender <- as.factor(nels$sex)
math.gender <- aggregate(math ~ homewrk + gender, data=nels, FUN="mean")

ggplot(math.gender, aes(x=homewrk, y=math, col=gender, type='l'),
            main='Regressions all in one plot')  + geom_line() 
text(0,77,'1=Male,',col='red')
text(1.5,77,'2=Female',col='lightseagreen')
```

But this works
```{r works_for_page36}
nels$gender <- as.factor(nels$sex)
math.gender <- aggregate(math ~ homewrk + gender, data=nels, FUN="mean")

plot(math.gender$homewrk,math.gender$math,type="n",
     main="Regression of Math on Homework \n Seperate lines for genders",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score")
lines(math.gender$homewrk[1:8],math.gender$math[1:8],col="blue")
lines(math.gender$homewrk[9:15],math.gender$math[9:15],col="red")
legend("topleft",legend=c("Female", "Male"),
       col=c("blue","red"), lty=1, cex=1.2)   
```

## Cross-classification with gender

This is around page 37
```{r cross_classification_gender}
with(nels,table(sex, homewrk))
```

## graph math ~ homework  + gender
These lines extend beyond data
```{r mean_homework_gender_beyond_data}
plot(math.gender$homewrk,math.gender$math,type="n",
     main="Regression of Mean Math on Homework \n Seperate lines for genders",
	 xlab="Time Spent Doing Homework",
	 ylab="Mean Math Score",
	 ylim=c(45,70))
abline(lm(math[1:7]~homewrk[1:7],data=nels),col="blue")
abline(lm(math[8:14]~homewrk[8:14],data=nels),col="red")
legend("topleft",legend=c("Female", "Male"),
       col=c("blue","red"), lty=1, cex=1.2)   
```

These do not extend beyond the range of data and fits lines to uncollapsed data
```{r just_data_math_homework_gender}
#  --- set up graphic frame:
plot(nels$homew,nels$math,
	type="n", 
	col="blue",
	lwd=2,
	main="NELS: Linear Regression by School by Gender \n(for where there is data)",
	xlab="Time Spent Doing Homework",
	ylab="Math Scores"
	)
# --- loop through the schools and add line
for (j in nels$sex) {
  sub <- nels[ which(nels$sex==j),]
  fitted <- fitted(lm(math~homew,sub))
  lines(sub$homew,fitted,col=j)
}  
```	

# SES


## page 41-42 Lattice plot with regressions
```{r page41_42}
xyplot(math ~ ses | school, data=nels, col.line='black', 
                      type=c('p','r'),
                      main='Varability in Math ~ ses relationship')
``` 
 
## page 43-44  Plot 3:  Lattice plot with smooth regressions  (i.e., loess)
```{r page43_44}
xyplot(math ~ ses | school, data=nels, col.line='black', 
                      type=c('p','smooth'),
                      main='Varability in Math ~ ses relationship')
```

## page 45  All regressions in one figure

```{r page45}
#  --- set up graphic frame:
par(ask=FALSE)
plot(nels$homew,nels$ses,
	type="n", 
	col="blue",
	lwd=2,
	main="NELS: Linear Regression by School \n(for where there is data)",
	xlab="Socio-econmic status (SES)",
	ylab="Math Scores",
	ylim=c(30,70),
	xlim=c(-2.5,2.5)
	)
# --- loop through the schools and add line
for (j in unique(nels$school)) {
  sub <- nels[ which(nels$school==j),]
  sub2 <- sub[order(sub$ses),]
  fitted <- fitted(lm(math~ses,sub2))
  lines(sub2$ses,fitted,col=j)
}     
```

## Create 10 groups with about the same number per group
```{r mean_mean_sesgroup}
nels$ses.grp <- as.numeric(cut(nels$ses, 10))
nels$ses.grp <- as.factor(nels$ses.grp)

mean.ses.grp <- aggregate(ses~ses.grp, data=nels, FUN="mean")
mean.math.grp<- aggregate(math~ses.grp, data=nels, FUN="mean")

# Plot grouped data
plot(mean.ses.grp$ses,mean.math.grp$math, type='b',
     main="Grouped SES:  math mean vs ses mean",
	 xlab="Mean SES",
	 ylab="Mean Math",
	 col="blue"	 )

# Plot math vs SES
plot(nels$ses,nels$math,type='p')	 
```

# Fitting models

First make categorical as factors.
```{r as.factors.are}
nels$white <- as.factor(nels$white)
nels$sex <- as.factor(nels$sex)
nels$public <- as.factor(nels$public)
```
## Empty/baseline

This is just a simple random intercept model.
```{r empty.model0}
summary(model0 <- lmer(math ~ (1  | school), data=nels, REML=FALSE,control = lmerControl(optimizer ="Nelder_Mead")))
(icc <- 24.85/(24.85+81.24))
bic.hlm(model0, nels$school)
```

## Preliminary lmm
The equations for this model can be found on pages 134--135 and the output on page 140.
```{r preliminaray.lmm}
summary( model1 <- lmer(math ~ homew + schCses + sex + white + schCses*white + schMses + public 
                 + (1 + homew | school), data=nels, REML=FALSE,
				 control = lmerControl(optimizer ="Nelder_Mead")))
screenreg(model1, single.row=TRUE)
bic.hlm(model1, nels$school)
```
## Too Complex

This model is too complex and does not converge.  It has an additional random slope.

### Random 1 + homew + schCses
```{r model1b}
model1b <- lmer(math ~ homew +  schCses + sex + white + schCses*white + schMses + public  + (1 + homew + schCses| school), data=nels, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
### Random 1 + homew + sex
```{r model1c}
model1c <- lmer(math ~ homew + schCses + sex + white + schCses*white + schMses + public + (1 + homew + sex| school), data=nels, REML=FALSE, control = lmerControl(optimizer="Nelder_Mead"))
```
### Random 1 + homew + white
This converges but has larger BIC.new and AIC compared to our prelininary model.
```{r model1d}
summary( model1d <- lmer(math ~ homew + schCses + sex + white + schCses*white + schMses + public  + (1 + homew + white| school), data=nels, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")))
```
## Test Random Effects (i.e., $\tau$s)	

### $H_o:  \tau_1^2 = \tau_{10}=0$

The full model is model 1 (the preliminary model) and we need to fit a null model for this test (i.e., remove ``homew'' from random part of model)
```{r model2}
summary(model2 <- lmer(math ~ homew + schCses + sex + white + schCses*white + schMses + public + (1 | school), data=nels, REML=FALSE,control = lmerControl(optimizer ="Nelder_Mead")))
```
The test statistic equals the difference between -2lnlike (i.e., deviances)
```{r mixture1}
(lr <- deviance(model2) - deviance(model1)  )
pvalue <- .5*(pchisq(lr,1,lower.tail=FALSE)) + .5*(pchisq(lr,2,lower.tail=FALSE))
round(pvalue, digits=2)
```

## Test fixed effects

### Mean SES vs sector
Since sector not significant in model 1, we'll try remove schMses (the predcitor of random intercept )
```{r model3}
summary(model3 <- lmer(math ~ homew + schCses + sex + white + schCses*white  + public + (1 + homew | school),        data=nels, REML=FALSE,
control =lmerControl(optimizer="Nelder_Mead")))

# for latex in lecture notes: texreg(list(model3), single.row=TRUE,custom.model.names("No Mean SES"), caption=("Why Secter was ns?"))
anova(model3, model1)
bic.hlm(model3, nels$school)
bic.hlm(model1, nels$school)
```

### For LR test for gender & public

```{r model1.reduced}
summary(model1.reduced <- lmer(math ~ homew + schCses + white + schCses*white + schMses + (1 + homew | school), data=nels, REML=FALSE,
control = lmerControl(optimizer="Nelder_Mead")))

anova(model1.reduced, model1)
bic.hlm(model1.reduced, nels$school)
# for lecture notes: texreg(list(model1.reduced),single.row=TRUE, custom.model.names=c("Simpler Model"))
deviance(model1.reduced)
```

### Drop sector
```{r drop.sector}
# page 156 
summary(model4 <- lmer(math ~ homew + schCses + sex + white + schCses*white  + schMses + (1 + homew | school),  data=nels, REML=FALSE,control = lmerControl(optimizer ="Nelder_Mead")))
```

### Drop interaction
Betwoen gender and white, which has not been significant
```{r no.interaction.model5}
summary(model5 <- lmer(math ~ homew + schCses + white + schMses + (1 + homew | school), data=nels, REML=FALSE,control = lmerControl(optimizer ="Nelder_Mead")))
```
And some checks
```{r post.model5}
anova(model5, model1.reduced)
# For laTex: texreg(list(model5),single.row=TRUE,custom.model.names=c("Final Model"))
bic.hlm(model5,nels$school)
```
## summary of model fit to data
```{r comparisons}
screenreg(
  list(model0, model1, model2, model3, model4, model5))
```
# Diagnositics for Model 5

Since model 5 at this points looks like it might be our final model, we'll do some diagnostic checks.

## Graph EBLUPS

Well do this in steps

### Step 1
We need to extract them from model object
```{r extract.Us}
ranU <- as.data.frame(ranef(model5))
```

### Step 2
Put them in separate objects
```{r create.u0j.u1j}
U0j <- ranU[which(ranU$term=="(Intercept)"), ]
U1j <- ranU[which(ranU$term=="homew"), ]
```

### Step 3
Sort rows by ranU, which is "condval"
```{r sort.condval}
U0j <-U0j[order(U0j$condval), ]
U1j <-U1j[order(U1j$condval), ]
```

### step 4 
Add integers for schools so that I can plot in order
```{r add.integer.id}
U0j$xgrp <- seq(1:23)
U1j$xgrp <- seq(1:23)
```
### Step 5 
Draw graph for U0j
```{r plot.U0j}
plot(U0j$xgrp, U0j$condval,type="p",
     main="Final Model: Random intercepts +/- 1 std error",
	 xlab="Schools (sort by Uo)",
	 ylab="EBLUP of U0j",
 col="blue"
	 )

# Add in the lines for sd errors by
# drawing arrows but with no arrowhead.
# This is a hack solution but it works
arrows(U0j$xgrp, U0j$condval-U0j$condsd,  
       U0j$xgrp, U0j$condval+U0j$condsd, 
	     length=0.00)
```
### Now for U1j:
Repeat above but for the random slopes
```{r graph.U1j}
plot(U1j$xgrp, U1j$condval, type="p",
   main="Final Model: Random slopes +/- 1 std error",
	 xlab="Schools (sort by U1j)",
	 ylab="EBLUP of U1j",
	 col="blue")

# For sd errors, draw arrows but with no arrowhead
arrows(U1j$xgrp, U1j$condval-U1j$condsd,  
       U1j$xgrp, U1j$condval+U1j$condsd, 
	     length=0.00)
```

## QQplot for U's

I will make custom qqplots that include +/- 1 sd error.
### Uoj
```{r qqplot.Uoj}
qqdata <- qqnorm(U0j$condval)

plot(qqdata$x, qqdata$y, type="p", pch=19, col="red",
main="Normal QQ plot of U0j with 95% Confidence Bands",
     xlab="Theoretical Value",
     ylab="Estimated U0j (EBLUPS)",
     ylim=c(-15,15))
arrows(qqdata$x, qqdata$y - 1.96*U0j$condsd, 
       qqdata$x, qqdata$y + 1.96*U0j$condsd, 
	     length=0.00, col="gray")
qqline(U0j$condval,col="steelblue")
```
### For U1j{
```{r qqplot.U1j}
qqdata <- qqnorm(U1j$condval)
plot(qqdata$x,qqdata$y,type="p",pch=19,col="red",
main="Normal QQ plot of U1j with 95% Confidence Bands",
     xlab="Theoretical Value",
     ylab="Estimated Uuj (EBLUPS)",
     ylim=c(-8,8)
     )
arrows(qqdata$x, qqdata$y - 1.96*U1j$condsd, 
       qqdata$x, qqdata$y + 1.96*U1j$condsd, 
	     length=0.00, col="gray")
qqline(U1j$condval,col="steelblue")
```
## Residual Plots

This gives standardized level 1 residuals.  These are stadnardized and 
```{r get.residuals}
res5 <- hlm_resid(model5, standardize=TRUE, include.ls=FALSE)
head(res5)
```

```{r plot.level1.conditional.residuals}
par(mfrow=c(2,2))
fit <- fitted(model5)
plot(fit, res5$.std.resid, 
     xlab='Conditional Fitted Values',
     ylab='Std Residuals',
     main='Conditional Level1 Residuals')
lines(c(min(fit),max(fit)),c(0,0), col="red", lwd=2)


# qqplot
qqnorm(res5$.std.resid)                    # draws plot
qqline(res5$.std.resid, col='red', lwd=2)  # reference line


# Histogram of standardized residuals with normal overlay
h<- hist(res5$.std.resid, breaks=15, density=20, col="seagreen")   # draw historgram
xfit <- seq(-3, 3, length=50)       # sets range & number of quantiles
yfit <- dnorm(xfit, mean=0, sd=1)   # should be normal(0,1)
yfit <- yfit*diff(h$mids[1:2])*length(res5$.std.resid) # use mid-points 
lines(xfit, yfit, col='blue', lwd=2)                   # draws normal


# Some things to describe the model
plot.new( ) 
text(.5,1.0,'Model 5')
text(.5,0.85,'Devience=3600.1')
text(.5,0.6,'AIC=3618.1')
text(.5,0.4,'BIC=3656.4')
```

## Level 2 or marginal
These are estimated random effects, which we got above using ranef().

```{r higher-level_residuals}
res5.margin <- resid_marginal(model5, type="pearson")
head(res5.margin)
```