# pacchetto
library(tidyverse)

# boxplot
npk %>% 
  ggplot(aes(block, yield)) + 
  geom_boxplot()

# ANOVA A DUE VIE
fit1 <- aov(yield ~ block * N, data = npk)

summary(fit1)

# oppure
library(broom)
(aov.df <- tidy(fit1))   # con la parentesi, salvo e visualizzo l'oggetto

# coefficienti contrasti
coef(fit1)

# tabella anova con contrasti
summary(fit1, 
        split = list(block = list(1,2,3,4,5)))

# ANOVA FATTORIALE
fit2 <- aov(yield ~ block + N * P * K, npk)

summary(fit2)

fit2