Plan <- read.table("https://husson.github.io/livre_exo/chap12_ex2.csv", header = TRUE, sep=";", encoding = "latin1") for (i in 1:3) Plan[,i] <- as.factor(Plan[,i]) summary(aov(Y ~ A + B + C, data=Plan)) X <- model.matrix(~ A+B+C, Plan, contrasts = "contr.sum") t(X)%*%X solve(t(X)%*%X)