This Rmd file shows the code to fit the models and produce the figures for the first part of the workshop. 

# Statistical model

## Model 1

```{r}
library(lme4)
library(tidyverse)
library(emmeans)
library(latex2exp)
```

```{r echo=TRUE}
url <- "https://raw.githubusercontent.com/k-state-id3a/mixed-models-fall25/refs/heads/newpart3/data/nitrogen_yield.csv"
df_n_ss  <-  read.csv(url)

m1 <- lm(Yield_SY ~ Total_N + I(Total_N^2), data = df_n_ss)
m2 <- lmer(Yield_SY ~ Total_N + I(Total_N^2) + (1|sy), data = df_n_ss)

m1
m2
```

```{r}
getME(m2, "X")
getME(m2, "Z")
getME(m2, "theta")
getME(m2, "b")
```

## Model 2 

```{r}
url <- "https://raw.githubusercontent.com/k-state-id3a/mixed-models-fall25/refs/heads/newpart3/data/cochrancox_kfert.csv"
df <- read.csv(url)

df$rep <- as.factor(df$rep)
df$K2O_lbac <- as.factor(df$K2O_lbac)
m_fixed <- lm(yield ~ K2O_lbac + rep, data = df)
m_random <- lmer(yield ~ K2O_lbac + (1|rep), data = df)

(mg_means_fixed <- emmeans(m_fixed, ~K2O_lbac, contr = list(c(1, 0, 0, 0, -1))))
(mg_means_random <- emmeans(m_random, ~K2O_lbac, contr = list(c(1, 0, 0, 0, -1))))
```


# Figures

```{r}
library(tidyverse)
library(MASS)
library(latex2exp)
```

## Figure 1 - attendees 

```{r }
d <- read.csv("mixedmodels.csv") %>% 
  mutate(Unit.or.department.at.K.State =
           str_replace(Unit.or.department.at.K.State, "Department of ", ""))

d <- d %>% 
  mutate(dept = case_when(Unit.or.department.at.K.State
                          %in% c("Animal Science", "Animal Science and Industry", 
                                 "Animal sciences and Industry", "ASI", "ASI ",
                                 "Animal Sciences and Industry", 
                                 "Animal science and industry ") ~ "Animal Science", 
                          Unit.or.department.at.K.State
                          %in% c("Department of Biological and Agricultural Engineering",
                                 "Biological and Agricultural Engineering",
                                 "Biological and Agricultural Engineering ",
                                 "biological.and agricultural engineering", 
                                 "Carl and Melinda dept of BAE", "BAE") ~ 
                            "Biological and Agricultural Engineering", 
                          Unit.or.department.at.K.State
                          %in% c("Plant Pathology ","Plant Pathology", "Plant pathology ", 
                                 "Plant pathology",
                                 "Plant Pathology/ Division of Biology ") ~ 
                            "Plant Pathology", 
                          Unit.or.department.at.K.State
                          %in% c("Agronomy","agronomy",
                                 "Wheat breeding/Agronomy", "Agronomy ") ~
                            "Agronomy", 
                          Unit.or.department.at.K.State
                          %in% c("K-state Libraries","K-state libraries", 
                                 "k-state libraries", "K-state libraries") ~
                            "K-State Libraries", 
                          Unit.or.department.at.K.State
                          %in% c("Biology","Biology ", "Division of Biology") ~ 
                            "Biology",
                          Unit.or.department.at.K.State
                          %in% c("Clinical sciences", "Clinical Sciences") ~ 
                            "Clinical Sciences",
                          Unit.or.department.at.K.State
                          %in% c("Center for Outcomes Research and Epidemiology,
                                 College of Veterinary Medicine") ~ "CORE", 
                          Unit.or.department.at.K.State
                          %in% c("Beef Cattle Institute, Diagnostic Medicine/Pathobiology",
                                 "Diagnostic Medicine/Pathobiology") ~
                            "Diagnostic Medicine/Pathobiology", 
                          TRUE ~ Unit.or.department.at.K.State) )

d %>%
  filter(dept != "") %>% 
  ggplot(aes(y = dept))+
  geom_bar(fill = "#B388EB")+
  theme_bw()+
  theme(panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        axis.title.y = element_blank())
# ggsave("mixed-models-fall25/images/day1/attendees.jpg")
```

## Figure 2 -- scatterplot 

```{r}
df_n_ss$yhat <- predict(m1)

df_fig <- expand.grid(
  Total_N = 0:max(df_n_ss$Total_N),
  sy = unique(df_n_ss$sy),
  Block = 1:4)
df_fig$yhat <- predict(m, newdata = df_fig)


df_n_ss %>% 
  ggplot(aes(Total_N, Yield_SY))+
  theme_bw()+
  geom_point(aes(), fill = "grey30",
             size = 3, alpha = .8,
             show.legend = F, shape = 21)+
  geom_line(aes(x = Total_N, y = yhat), data = df_fig)+
  labs(x = expression(N~fertilizer~(lb~N~ac^{-1})), 
       y = expression(Yield~(kg~ha^{-1})))+
  theme_classic()
# ggsave("mixed-models-fall25/images/day1/linearmodel1.jpg", width = 7, height = 5)
```

## Figure 3 -- scatterplot quadratic model 2 (with location-specific intercepts)

```{r}
df_fig$yhat.mixed <- predict(m2, newdata = df_fig)

df_n_ss %>% 
  ggplot(aes(Total_N, Yield_SY))+
  theme_bw()+
  rcartocolor::scale_fill_carto_d(palette = "Vivid")+
  rcartocolor::scale_color_carto_d(palette = "Vivid")+
  geom_point(aes(fill = factor(sy)),
             size = 3, alpha = .8,
             show.legend = F, shape = 21)+
  geom_line(aes(x = Total_N, y = yhat.mixed, 
                group = paste(sy, Block),
                color = factor(sy)), 
            show.legend = F,
            data = df_fig)+
  labs(x = expression(N~fertilizer~(lb~N~ac^{-1})), 
       y = expression(Yield~(kg~ha^{-1})))+
  theme_classic()
# ggsave("mixed-models-fall25/images/day1/linearmodel2.jpg", width = 7, height = 5)

```


## Figure 4 -- Normal distributions 


```{r echo=FALSE, fig.width=10, message=FALSE, warning=FALSE, fig.align='center'}
set.seed(42)

x <- seq(-6, 6, length=100)
hx <- dnorm(x)
labels <- c("x ~ Normal(0, 1)", 
            "x ~ Normal(0, 16)", "x ~ Normal(2, 1)")
colors <- c("darkgreen", "tomato", "purple")


# png(filename="../../images/day1/normal_univariate.png")
plot(x, dnorm(x, 0, 1), lwd=2, col= "darkgreen",
     type="l", lty= 1, xlab="x",
  ylab="[x]", main="Univariate Normal")

# lines(x, dnorm(x, 0, 1), lwd=2, col= "darkgreen")
lines(x, dnorm(x, 0, 4), lty = 2, lwd=2, col= "tomato")
lines(x, dnorm(x, 2, 1), lty = 3, lwd=2, col= "purple")

legend("topleft", inset=.05, title="Distributions",
  labels, lwd=2, lty=c(1, 2, 3), col=colors)
# dev.off()
```

## Figure 5 -- Multivariate normal distribution  

```{r}
set.seed(42)
dat <- mgcv::rmvn(1000, c(10, 8), matrix(c(1, .6, .6, 1), 2, 2))

png(filename="../../images/day1/normal_multivariate.jpg")
plot(dat[,1], dat[,2], xlab = TeX("$y_1$"), ylab = TeX("$y_2$"), main = "")
segments(x0 = 10, y0 = 4.5, x1 = 10, y1 = 8, col = "tomato", lwd = 3)
segments(x0 = 6.5, y0 = 8, x1 = 10, y1 = 8, col = "darkgreen", lwd = 3)
lines(sort(dat[,1]), 5 + 5*dnorm(sort(dat[,1]), 10, 1), lty = 2, lwd=3, 
      col= "tomato")
lines(7 + 2 * dnorm(sort(dat[which(dat[,2]<8),2]), 8, 1), 
      sort(dat[which(dat[,2]<8),2]),
      lty = 3, lwd=3, col= "darkgreen")
lines(7 + 2 * dnorm(sort(dat[which(dat[,2]>8),2]), 8, 1), 
      sort(dat[which(dat[,2]>8),2]),
      lty = 3, lwd=3, col= "darkgreen")
dev.off()
```

## Figure from the button -- Multivariate normal distribution with zero covariance  

```{r}
set.seed(42)
dat <- mgcv::rmvn(1000, c(10, 8), matrix(c(1, 0, 0, 1), 2, 2))

png(filename="../images/day1/normal_multivariate_0cov.jpg")
plot(dat[,1], dat[,2], xlab = TeX("$y_1$"), ylab = TeX("$y_2$"), main = "")
segments(x0 = 10, y0 = 4.5, x1 = 10, y1 = 8, col = "tomato", lwd = 3)
segments(x0 = 6.5, y0 = 8, x1 = 10, y1 = 8, col = "darkgreen", lwd = 3)
lines(sort(dat[,1]), 5 + 5*dnorm(sort(dat[,1]), 10, 1), lty = 2, lwd=3, 
      col= "tomato")
lines(7 + 2 * dnorm(sort(dat[which(dat[,2]<8),2]), 8, 1), 
      sort(dat[which(dat[,2]<8),2]),
      lty = 3, lwd=3, col= "darkgreen")
lines(7 + 2 * dnorm(sort(dat[which(dat[,2]>8),2]), 8, 1), 
      sort(dat[which(dat[,2]>8),2]),
      lty = 3, lwd=3, col= "darkgreen")
dev.off()
```

## Figure 6 -- Visual representation of the variance-covariance matrix assuming independent observations 

```{r}
diag(sigma(m)^2, nrow = nrow(df_n_ss)) %>% 
  as.data.frame() %>% 
  rownames_to_column("x") %>% 
  pivot_longer(cols = V1:V158, names_to = "y") %>%
  mutate(y = str_replace(y, "V", "") %>% as.numeric(),
         x = as.numeric(x)) %>% 
  ggplot(aes(x,-y))+
  geom_tile(aes(fill = value))+
  theme_minimal()+
  scale_fill_viridis_c()+
  theme(aspect.ratio = 1, 
        axis.text = element_blank(), 
        legend.position = "none", 
        panel.grid = element_blank(), 
        axis.title = element_blank(), 
        axis.ticks = element_blank())
# ggsave("mixed-models-fall25/images/day1/covmatrix_crd.jpg", width = 7, height = 5)
```

## Figure 7 -- Visual representation of the variance-covariance matrix from a mixed model

```{r}
Z <-  getME(m2, "Z") 
sigma2_s <- as.data.frame(VarCorr(m2))[1,4]
G <- diag(c(rep(sigma2_s, 5)))

Z %*% G %*% t(Z) %>%
  as.matrix() %>% 
  as.data.frame() %>% 
  rownames_to_column("x") %>% 
  pivot_longer(cols = `1`:`158`, names_to = "y") %>%
  mutate(x = as.numeric(x), y= as.numeric(y)) %>% 
  
  left_join(
  diag(sigma(m)^2, nrow = nrow(df_n_ss)) %>% 
  as.data.frame() %>% 
  rownames_to_column("x") %>% 
  pivot_longer(cols = V1:V158, names_to = "y") %>%
  mutate(y = str_replace(y, "V", "") %>% as.numeric(),
         x = as.numeric(x)) %>% 
  rename(s2 = value)) %>% 

  ggplot(aes(-x,y))+
  geom_tile(aes(fill = value+s2))+
  theme_minimal()+
  scale_fill_viridis_c()+
  theme(aspect.ratio = 1, 
        axis.text = element_blank(), 
        legend.position = "none", 
        panel.grid = element_blank(), 
        axis.title = element_blank(), 
        axis.ticks = element_blank())
# ggsave("mixed-models-fall25/images/day1/covmatrix_rcbd.jpg", width = 7, height = 5)
```

## Figure 8 -- RCBD blocks fixed/random


```{r}
as.data.frame(mg_means_fixed$emmeans) %>%
  mutate(blocks = "fixed") %>% 
  bind_rows(as.data.frame(mg_means_random$emmeans) %>%
              mutate(blocks = "random")) %>% 
  ggplot(aes(K2O_lbac, emmean))+
  geom_errorbar(aes(ymin = emmean-SE, ymax = emmean+SE, 
                    color = blocks), 
                position = position_dodge(width = .2),
                width = 0)+
  geom_text(aes(x = 4.6, y = 8), label = "s.e.(difference\nbetween means)",
            color = "grey30", size = 3.5)+ 
  geom_point(aes(x = 5.2, y = 8), color = "grey50")+ 
  geom_point(aes(x = 5.4, y = 8), color = "grey50")+ 
  geom_errorbar(aes(x = 5.4, y = 8, ymin = 8 - SE, ymax = 8 + SE), 
                width = 0, color = "grey50",
                data = as.data.frame(mg_means_random$contrasts))+
  geom_errorbar(aes(x = 5.2, y = 8, ymin = 8 - SE, ymax = 8 + SE), 
                width = 0, color = "grey50", 
                data = as.data.frame(mg_means_fixed$contrasts))+
  scale_color_manual(values = c("grey30", "tomato"))+
  
  geom_point(aes(color = blocks), position = position_dodge(width = .2))+
  theme_classic()+
  labs(title = "Spoiler: difference in results for RCBD data\nmodeled with fixed vs. random blocks", 
       color = "Blocks modeled as",
       y = expression(Yield~(tn~ac^{-1})),
       x = expression(K[2]~O~(lb~ac^{-1})))+
  theme(legend.position = "bottom", 
        plot.title = element_text(hjust = .5))
# ggsave("mixed-models-fall25/images/day1/emmeans_blocks.jpg", width = 6, height = 4)
```


## Figure 9 -- shrinkage: blocks fixed/random

```{r}
data.frame(b0_j = c(as.vector(getME(m2, "b")), as.data.frame(emmeans(m2_fixed, ~sy, at = list(Total_N = 0)))$emmean - mean(as.data.frame(emmeans(m2_fixed, ~sy, at = list(Total_N = 0)))$emmean)), 
           model = rep(c("random", "fixed"), each = 5)) %>% 
  ggplot(aes(x = b0_j , y = 0))+
  geom_vline(xintercept = 0, linetype = 2, color = "tomato")+
  stat_function(fun = dnorm, 
                args = list(mean = 0, sd = sqrt(sigma2_s)),
                color = "tomato", 
                xlim = c(-1750, 1750))+
  theme_classic()+
  scale_fill_manual(values = c("grey30", "tomato"))+
  labs(x = TeX("$\\hat{u}_j$ (Block effect)", italic =T), 
       y = TeX("$\\[ \\hat{u}_j \\]$", italic =T),
       fill = "Blocks modeled as")+
  theme(aspect.ratio = .25, legend.position = "bottom")+
  geom_point(bins = 20, aes(fill = model), shape=25, size =3, alpha = .8)+
  coord_cartesian(xlim = c(-1750, 1750))

# ggsave("mixed-models-fall25/images/day1/shrinkage.jpg", width = 6, height = 2.5)
```