calculating ANOVA 'by hand' I've calculated an ANOVA in R using the aov() function, then calculated it semi-manually using my own code. I'm interested in why my mean sum of squares between groups (0.634443) is not the same as that outputted by aov() (0.4145), and yet my mean sum of squares within groups is the same as that outputted by aov()?
# create fake data
set.seed(50)
x = data.frame(length=rnorm(9), 
               site = c(rep('a', 3), 
                        rep('b', 3), 
                        rep('c', 3)), stringsAsFactors=F)
# ANOVA summary
summary(aov(length ~ factor(site), x))
# Df Sum Sq Mean Sq F value Pr(>F)
# factor(site)  2  0.829  0.4145   0.514  0.622
# Residuals     6  4.840  0.8067 

# calculate sum of squares within, between and total
sumlength <- sum(x$length)   
    n = 3
    overallMean <- sumlength/n
    aMean <- mean(x$length[x$site=="a"])
    bMean <- mean(x$length[x$site=="b"])
    cMean <- mean(x$length[x$site=="c"]) 

ssw <- sum((x$length[x$site == 'a'] - aMean) * (x$length[x$site == 'a'] - aMean)) + 
            sum((x$length[x$site == 'b'] - bMean) * (x$length[x$site == 'b'] - 
            bMean)) + sum((x$length[x$site == 'c'] - cMean) * (x$length[x$site == 
            'c'] - cMean))  
sst <- sum((x$length - overallMean) * (x$length - overallMean))  
ssb <- ((aMean - overallMean) * (aMean - overallMean) + (bMean - overallMean) * 
        (bMean - overallMean) + (cMean - overallMean) * (cMean - overallMean)) * 3

meanSsb <- ssb/2     # 0.634443
meanSsw <- ssw/6     # 0.8066979

 A: output:
summary(aov(length ~ factor(site), x))
             Df Sum Sq Mean Sq F value Pr(>F)
factor(site)  2  0.829  0.4145   0.514  0.622
Residuals     6  4.840  0.8067    


overallmean should be the overall mean

overallmean <- mean(x$length)

(you divided by 3 rather than dividing by 9).

To get the sum of squares for site, you need to compute
$$
\sum_i n_i \,(\bar{Y}_{i\cdot} - \bar{Y}_{\cdot\cdot})^2
$$
i.e., the sum of weighted (with weights equal to $n_1, n_2, \dotsc$) squared deviations of each estimated factor level mean ($\bar{Y}_{i\cdot}$) around the overall mean ($\bar{Y}_{\cdot\cdot}$). 
In R, (using the same kind of calculations as you)
na <- nrow(x[x$site == "a", ])
nb <- nrow(x[x$site == "b", ])
nc <- nrow(x[x$site == "c", ])

aMean <- mean(x[x$site == "a", "length"])
bMean <- mean(x[x$site == "b", "length"])
cMean <- mean(x[x$site == "c", "length"])

sum(na * (aMean - overallmean)^2 + 
    nb * (bMean - overallmean)^2 + 
    nc * (cMean - overallmean)^2)
[1] 0.8290974


To get the residual sum of squares, you need 
$$
\sum_i \sum_j (Y_{ij} - \bar{Y}_{i\cdot})^2
$$
i.e., the sum of squared deviations of each observation ($Y_{ij}$) around its estimated factor level mean ($\bar{Y}_{i\cdot}$)
sum(c((x[x$site == "a", "length"] - aMean)^2,
      (x[x$site == "b", "length"] - bMean)^2,
      (x[x$site == "c", "length"] - cMean)^2))
[1] 4.840187

A: Look at the line n = 3.  The overall number of data points is 9, not 3, so when you calculate the overall mean you are dividing the sum of 9 numbers by 3, so your mean is not the mean.  Change that to n = 9 and it looks like it matches to me.
