---
title: "STAT 516 Lec 05 in class"
format: html
---

```{r}
link <- url("https://people.stat.sc.edu/gregorkb/data/KNNLrust.txt")
rust <- read.csv(link,col.names=c("score","brand","rep"),sep = "", header = FALSE)
head(rust)
```


```{r}
boxplot(score ~ as.factor(brand), data = rust)
```

```{r}
lm_out <- lm(score ~ as.factor(brand), data = rust)
lm_out
```

Compute all the quantities in the ANOVA table:

```{r}
Y <- rust$score
Y.. <- mean(Y) # overall mean
Yi. <- aggregate(score ~ brand,mean, data = rust)[,2]
n <- 10
a <- 4
N <- a*n

SStot <- sum( (Y - Y..)**2)
SStrt <- sum( n*(Yi. - Y..)^2)
SSerror <- SStot - SStrt

MStrt <- SStrt / (a - 1)
MSerror <- SSerror / (N - a)

Ftest <- MStrt / MSerror
Ftest
alpha <- 0.05
Fcrit <- qf(1 - alpha, a - 1, N - a)
Fcrit

pval <- 1 - pf(Ftest, a - 1, N - a)
```

An easier way:

```{r}
lm_out <- lm(score ~ as.factor(brand), data = rust)
anova(lm_out)
```

```{r}
aggregate(score ~ brand,mean, data = rust)
```

Tukey interval using value from the table:
```{r}
lo <- 43.14 - 89.44 - 3.85 * sqrt(6.1)/sqrt(10)
up <- 43.14 - 89.44 + 3.85 * sqrt(6.1)/sqrt(10)
c(lo,up)
```

Tukey interval using exact value:
```{r}
lo <- 43.14 - 89.44 - qtukey(1-alpha,a,a*(n-1)) * sqrt(6.1)/sqrt(10)
up <- 43.14 - 89.44 + qtukey(1-alpha,a,a*(n-1))* sqrt(6.1)/sqrt(10)
c(lo,up)
```



Tukey's automatically:

```{r}
aov_out <- aov(score ~ as.factor(brand), data = rust)
TukeyHSD(aov_out)
plot(TukeyHSD(aov_out))
```


Dunnett's with group 1 as baseline with approximate number from table: Group 2 minus group 1.

```{r}
lo <- 89.44 - 43.14 - 2.47 * sqrt(6.1)*sqrt(2/10)
up <- 89.44 - 43.14 + 2.47 * sqrt(6.1)*sqrt(2/10)
c(lo,up)
```

Dunnett's automatically:

```{r}
library(DescTools) # first time run install.packages("DescTool")
Dunnetts_out <- DunnettTest(score ~ as.factor(brand), data = rust, control = "1")
Dunnetts_out
plot(Dunnetts_out)
```

Check assumptions:
```{r}
plot(lm_out,which=1,add.smooth= F)

plot(lm_out, which = 2)
```

Levene's test for equal variances.

Get residuals from original one-way ANOVA model:
```{r}
eabs <- abs(lm_out$residuals)
lm_out2 <- lm(eabs ~ as.factor(rust$brand))

summary(lm_out2)
```

```{r}
library(car)
leveneTest(score ~ as.factor(brand), data = rust, center = mean)
```