---
title: "ANOVA and Such"
author: "Carolyn J. Anderson"
date: "12/28/2022"
output: 
  html_document:
    toc: yes
---

This example using NELS data but with 10 schools.  

# Set Up

First libraries that we'll use
```{r}
library(lme4)
library(lattice)
library(stargazer)
```

Read in data and take a look
```{r}
setwd("D:/Dropbox/edps587/lectures/2 anovas")

nels10 <- read.table(file="NELS10.txt", header=TRUE)
names(nels10)
head(nels10)
tail(nels10)
```
I tend to be a bit lazy, so can just use variable name without also using data table name
```{r}
attach(nels10)
```

# A First Look

## Descriptive Statistics
We could use base R
```{r}
summary(math)
mean(math)
var(math)

# Summary statistics by school

mn <- aggregate(math~sch_id, data=nels10,  mean)
sd <- aggregate(math~sch_id, data=nels10,  sd)
med <- aggregate(math~sch_id, data=nels10, median)
simple <- cbind(mn, sd[,2], med[,2])
round(simple,digits=2)
```

Or the stargazer package
```{r}
stargazer(nels10, type="text")
```
and by school,
```{r eval=FALSE}
by(nels10, nels10$sch_id, stargazer, type = 'text')
```

## Plot the data

We'll plot school mean maths scores with +/- one standard deviation.
```{r}
# -- a numbers 1 to 10 as index variable
simple$index <- seq(1:10)

plot(simple$index, simple$math, pch=19,
     xlab="School",
	 ylab="Math",
	 xlim=c(min(simple$index),max(simple$index)),
	 ylim=c(30,70),
	 main="School Mean Math Score +/- one Standard Deviation"
	 )
lines(rbind(1:10,1:10,NA),rbind(simple$math-simple[,3],simple$math+simple[,3],NA))
```

# Some Models

## No Predictors -- just the mean of math
This is just plain old linear regression (no predictors), which is just the overall mean for math.
```{r}
summary(linreg <- lm(math ~ 1, data=nels10))
anova(linreg)                    # to get ANOVA table
plot(linreg, las=1)              # plots fitted by residuals
```

## ANOVA -- least squares

Here the predictors are schools so that each school will have different mean

```{r}
id <- as.factor(sch_id)
summary( fixedANOVA <- lm(math ~ id, data=nels10,REML=TRUE) )
anova(fixedANOVA)
```

Note that the estimation method is REML, which will give same results as standard ANOVA

## ANOVA via MLE

Check to see if these are different results than those when using REML=TRUE.

```{r}
id <- as.factor(sch_id)
summary( fixedANOVA2 <- lm(math ~ id, data=nels10, REML=FALSE) )
anova(fixedANOVA2)
```

## Random Intercept HLM (Random Effects ANOVA))

REML is the default when using lmer

```{r}
summary( ranintREML <- lmer(math ~ 1 + (1 | sch_id), data=nels10, REML=TRUE) )
```

And when using MLE,
```{r}
summary( ranintMLE <- lmer(math ~ 1 + (1 | sch_id), data=nels10, REML=FALSE) )
```

# Intra-class Correlation

Just cutting and pasting from random intercept model
```{r}
(icc <- (30.54)/(30.54 + 72.24))
```

Or without cutting and pasting,
```{r icc}
s <- summary(ranintMLE)
# -- see what available...we want what's in "varcor" as a data frame
names(s)
(s <- as.data.frame(s$varcor))
(icc2 <- s$vcov[1] / (s$vcov[1] + s$vcov[2]))
```
Later we'll write a little function to do this.

# Linear Regressions 

First a model with an intercept and homework
```{r}
multreg <- lm(math ~ 1 + homework, data=nels10)
summary(multreg)
anova(multreg)
fitmath <- fitted(multreg)
plot(homework, math, 
      type='p', 
      main='Linear Regression of Math on Time Spent Doing Homework',
      xlab='Time Spent Doing Homework', 
      ylab='Fitted Math Scores from Linear Regression')
lines(homework,fitmath,col="blue")
```

# By School

But we should look at things for each school
```{r}
# Lattice plot with linear regressions
xyplot(math ~ homework | id, col.line='black', 
       type=c('p','r'),
       cex.xlab=1.5,
	   cex.ylab=1.5,
	   xlab="Time Spent Doing Homework",
	   ylab="Math Scores",
	   cex.main=1.2,
       main='Variability in Math ~ Homework relationship')
```

This will basically fit a linear regression to each school

```{r}
mathbysch <- lm(formula = math ~ 1*id + homework * id, data = nels10)
summary(mathbysch)
fschmath <- fitted(mathbysch)
```

We can plot results
```{r}
school.list <- unique(id)
length(school.list)

# Pull out intercepts and slopes from model objects
coef <- coefficients(mathbysch)
intercepts <- c(0, coef[3:11])
slopes <- c(0, coef[12:20])

# Plot regressions
plot(data = nels10, math~homework,type = 'n', 
     ylim = c(min(math),max(math)),
     xlim = c(min(homework),max(homework)),
     cex.main = 1.25,
     xlab = 'Time Spent Doing Homework', ylab = "Math Scores",
     main = "Separate Regression for Each School")
colors <- c("red", "orange", "darkgoldenrod1","magenta",
                      "green", "blue", "purple",
                      "white", "grey", "black")
for(i in 1:length(school.list)){
  abline(a = (coef[1] + intercepts[i]), b = (coef[2]+slopes[i]),
              col = colors[i])
  par<- par(new=F)
}
```

We'll use lmer and fit was is analogous to fitting regressions to each school,
but the school effects are random.   line is the average school and the multicolor 
are the conditional regressions (i.e., they include school specific estimates of 
intercept & slope)

```{r}
summary( model2 <- lmer(math ~ homework + ( homework | id), data=nels10, REML=FALSE) )

school.list <- unique(id)
length(school.list)

# Pull the fixed and random effects
fixed.effects <- fixef(model2)    # these are gammas
rand.effects <- ranef(model2)     # these are Us for schools
# these are used to get conditional regressions (conditional on school)
params.student <- cbind( (rand.effects$id[1]+fixed.effects[1]),
                         (rand.effects$id[2]+fixed.effects[2]) )

plot(data = nels10, math~homework,type = 'n', 
     ylim = c(min(math),max(math)),
     xlim = c(min(homework),max(homework)),
     cex.main = .75,
     xlab = 'Time Spent Doing Homework', ylab = "Math Scores",
     main = "Estimated Conditional Models from Random Intercepts & Slope Model")

abline(a=fixed.effects[1], b= fixed.effects[2], col='black', lwd=2)

for(i in 1:length(school.list)){
  abline(a = params.student[i,1], b = params.student[i,2],col = i)
  par<- par(new=F)
}
```