```
x.doublebar <- mean(c(a, b, c))
x.bar <- apply(cbind(a,b,c), 2, mean)
n.i <- apply(cbind(a,b,c), 2, length)
ssd.between <- sum(n.i*(x.bar-x.doublebar)^2)
ssd.within <- sum((rbind(a,b,c)-x.bar)^2)
ssd.total <- sum((c(a,b,c)-x.doublebar)^2)
df.between <- 3-1
df.within <- length(c(a,b,c))-3
ms.between <- ssd.between / df.between
ms.within <- ssd.within / df.within
F <- ms.between / ms.within
P <- 1-pf(F, df.between, df.within)
```

We get an F statistic of 51.6547 with 2 and 87 degrees of freedom and a P-value of 1.5543 × 10^{-15}.