```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Load R packages 
```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(glmmTMB)
library(agridat)
```

Fit models
```{r}
# load data
data("omer.sorghum")
df <- omer.sorghum |> filter(env == "E3")

# find out how many levels for trt factor "genotype"
n_distinct(df$gen)
# find out how many blocks
n_distinct(df$rep)

# fit linear model
m <- lm(yield ~ 1 + gen + rep, data= df)

# inspect X
model.matrix(m)

coef(m)

nrow(df) - 3 -17 -1
# these options will affect the baseline for the estimation of the effects 
# options(contrasts = c("contr.sum", "contr.poly"))
# options(contrasts = c("contr.treatment", "contr.poly"))
```

```{r}

options(contrasts = c("contr.sum", "contr.poly"))

m3 <- lm(yield ~ 1 + gen + rep, data= df)

# inspect X
model.matrix(m3)

coef(m3)

nrow(df) - 3 -17 -1
# these options will affect the baseline for the estimation of the effects 
# 
# options(contrasts = c("contr.treatment", "contr.poly"))
```



```{r}
# check for constant variance 
plot(m$fitted.values, m$residuals)

plot(m$fitted.values, abs(m$residuals))

# summary(lm(abs(m$residuals) ~ m$fitted.values))
# car::leveneTest(yield ~ gen , data = df)

# check for normality 
car::qqPlot(m)

# shapiro.test(m$residuals)

```

```{r}
car::Anova(m, type = 2)
```

```{r}
aov(yield ~ 1 + gen + Error(rep), data= df)
```

```{r}
#get sigma hat
sqrt(1306996/51)

sigma(m)
```


```{r}
emmeans::emmeans(m, ~gen)
```


```{r}
m2 <- lm(yield ~ 1 + gen, data= df)
car::Anova(m2, type = 2)
emmeans::emmeans(m2, ~gen)

sigma(m)^2;sigma(m2)^2
```


```{r}
df2 <- df |> filter(rep %in% c("R1", "R4"))

# find out how many levels for trt factor "genotype"
n_distinct(df2$gen)
# find out how many blocks
n_distinct(df2$rep)

# fit the same model with the new (smaller) data
m2 <- lm(yield ~ gen + rep, data= df2)

# check assumptions
plot(m2$fitted.values, abs(m2$residuals))
car::qqPlot(m2)

# see how DF of the error change!
car::Anova(m, type = 2)
car::Anova(m2, type = 2)
```