RCBD PRACTICE

detailed calculation

obs <- c(9, 11, 3, 7, 8, 13, 5, 10, 7, 12, 8, 4)
treat <- rep(c("A", "B", "C", "D"), times = 3)
block <- rep(c(1, 2, 3), each = 4)
data <- data.frame(obs, treat, block)
data

   obs treat block
1    9     A     1
2   11     B     1
3    3     C     1
4    7     D     1
5    8     A     2
6   13     B     2
7    5     C     2
8   10     D     2
9    7     A     3
10  12     B     3
11   8     C     3
12   4     D     3

a <- 4
b <- 3
N <- a*b

CF <- sum(obs)^2/N
CF

[1] 784.0833

SST <- sum(obs^2) - CF
SST

[1] 106.9167

treat_total <- tapply(obs, treat, sum)
treat_total

 A  B  C  D 
24 36 16 21 

SStreat <- sum(treat_total^2/b) - CF
SStreat

[1] 72.25

block_total <- tapply(obs, block, sum)
block_total

 1  2  3 
30 36 31 

SSblock <- sum(block_total^2/a) - CF
SSblock

[1] 5.166667

SSE <- SST - SStreat - SSblock
SSE

[1] 29.5

df1 <- a-1
df2 <- b-1
dfe <- df1 * df2
dft <- N-1

MSTreat <- SStreat/df1
MSTreat

[1] 24.08333

MSBlock <- SSblock/df2
MSBlock

[1] 2.583333

MSE <- SSE / dfe
MSE

[1] 4.916667

ftreat <- MSTreat/MSE
ftreat

[1] 4.898305

fblock <- MSBlock/MSE
fblock

[1] 0.5254237

ftreatab <- qf(0.98, df1, dfe, FALSE )
ftreatab

[1] 7.28698

fblocktab <- qf(0.98, df2, dfe, FALSE)
fblocktab

[1] 8.052094

if(ftreat > ftreatab){ 
  print("H0 may be rejected for treatment")
}else{
    print("H0 may not be rejected for treatment")
  }

[1] "H0 may not be rejected for treatment"

if(fblock > fblocktab){ 
  print("H0 may be rejected for blokcs")
}else{
    print("H0 may not be rejected for blokcs")
  }

[1] "H0 may not be rejected for blokcs"

shortcut

data$treat <- as.factor(data$treat)
data$block <- as.factor(data$block)

summary(aov(obs~treat))

            Df Sum Sq Mean Sq F value Pr(>F)  
treat        3  72.25  24.083   5.558 0.0234 *
Residuals    8  34.67   4.333                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(obs~block))

            Df Sum Sq Mean Sq F value Pr(>F)
block        1   0.12   0.125   0.012  0.916
Residuals   10 106.79  10.679               

analysis <- aov(obs~treat+block, data)
summary(analysis)

            Df Sum Sq Mean Sq F value Pr(>F)  
treat        3  72.25  24.083   4.898 0.0471 *
block        2   5.17   2.583   0.525 0.6162  
Residuals    6  29.50   4.917                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(analysis)

overall_mean <- mean(obs)
overall_mean

[1] 8.083333

treatmean <- tapply(obs, treat, mean)
treatmean

        A         B         C         D 
 8.000000 12.000000  5.333333  7.000000 

blockmean <- tapply(obs, block, mean)
blockmean

   1    2    3 
7.50 9.00 7.75 

blockeffect <- blockmean - overall_mean
blockeffect

         1          2          3 
-0.5833333  0.9166667 -0.3333333 

treateffect <- treatmean - overall_mean
treateffect

          A           B           C           D 
-0.08333333  3.91666667 -2.75000000 -1.08333333 

residuals(analysis)

          1           2           3           4           5           6 
 1.58333333 -0.41666667 -1.75000000  0.58333333 -0.91666667  0.08333333 
          7           8           9          10          11          12 
-1.25000000  2.08333333 -0.66666667  0.33333333  3.00000000 -2.66666667