CM2

# Simulate or load the Latin square dataset
# Example dataset
data <- data.frame(
  Treatment = factor(rep(c("A", "B", "C"), each = 3)),
  Block1 = factor(rep(c("Block1_1", "Block1_2", "Block1_3"), 3)),
  Block2 = factor(c("Block2_1", "Block2_2", "Block2_3", 
                    "Block2_2", "Block2_3", "Block2_1",
                    "Block2_3", "Block2_1", "Block2_2")),
  Response = c(15, 18, 20, 14, 19, 22, 13, 17, 21)
)
head(data)

##   Treatment   Block1   Block2 Response
## 1         A Block1_1 Block2_1       15
## 2         A Block1_2 Block2_2       18
## 3         A Block1_3 Block2_3       20
## 4         B Block1_1 Block2_2       14
## 5         B Block1_2 Block2_3       19
## 6         B Block1_3 Block2_1       22

# Load library
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Fit the model
latin_square_model <- aov(Response ~ Treatment + Block1 + Block2, data = data)

# Display ANOVA table
summary(latin_square_model)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## Treatment    2   2.67    1.33    1.00 0.5000  
## Block1       2  74.00   37.00   27.75 0.0348 *
## Block2       2   0.67    0.33    0.25 0.8000  
## Residuals    2   2.67    1.33                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Sum of squares calculations
SST <- sum((data$Response - mean(data$Response))^2)
SSTr <- sum((tapply(data$Response, data$Treatment, mean) - mean(data$Response))^2)
SSB1 <- sum((tapply(data$Response, data$Block1, mean) - mean(data$Response))^2)
SSB2 <- sum((tapply(data$Response, data$Block2, mean) - mean(data$Response))^2)
SSE <- SST - SSTr - SSB1 - SSB2

list(SST = SST, SSTr = SSTr, SSB1 = SSB1, SSB2 = SSB2, SSE = SSE)

## $SST
## [1] 80
## 
## $SSTr
## [1] 0.8888889
## 
## $SSB1
## [1] 24.66667
## 
## $SSB2
## [1] 0.2222222
## 
## $SSE
## [1] 54.22222

# F-ratio calculations
MS_T <- SSTr / (length(unique(data$Treatment)) - 1)
MS_B1 <- SSB1 / (length(unique(data$Block1)) - 1)
MS_B2 <- SSB2 / (length(unique(data$Block2)) - 1)
MSE <- SSE / ((length(unique(data$Treatment)) - 1) * 
              (length(unique(data$Treatment)) - 2))

F_T <- MS_T / MSE
F_B1 <- MS_B1 / MSE
F_B2 <- MS_B2 / MSE

list(F_T = F_T, F_B1 = F_B1, F_B2 = F_B2)

## $F_T
## [1] 0.01639344
## 
## $F_B1
## [1] 0.454918
## 
## $F_B2
## [1] 0.004098361