start.time = Sys.time()
set.seed(5)
no = 10
np = 20
nr = 4
d = 4
X = matrix(rnorm(no*np*nr*d),ncol = d)
X = data.frame(operator = rep(1:no,np*nr), part = rep(1:np, no,, nr), replication = rep(1:nr,,, no*np), X)
Mat_exp = expand.grid(1:no, 1:np)

Reduce('+', lapply(1:nrow(Mat_exp), function(i){
    M = t(as.matrix(X[,4:ncol(X)]))
    X_demean = M - rowMeans(M[,X$operator == Mat_exp[i,1]]) - rowMeans(M[,X$part == Mat_exp[i,2]]) + mean(M)
    X_demean %*% t(X_demean)
})) / (nrow(X) - no - np + nr)
Sys.time() - start.time

library(plyr)
library(data.table)
rm(list=ls(all=TRUE))
start.time = Sys.time()
set.seed(5)
no = 10
np = 20
nr = 4
d = 4
X = matrix(rnorm(no*np*nr*d),ncol = d)
X = data.table(operator = rep(1:no,np*nr), part = rep(1:np, no,, nr), replication = rep(1:nr,,, no*np), X)
setkey(X, operator,part)
Mat_exp = expand.grid(1:no, 1:np)
Reduce('+', alply(matrix(1:nrow(Mat_exp),,1),1, function(i){
    X_demean = apply(X[, 4:ncol(X), with = FALSE], 1, function(v) v-colMeans(X[J(1:no,Mat_exp[i,2]), 4:ncol(X), with = FALSE])-colMeans(X[J(Mat_exp[i,1],1:np), 4:ncol(X), with = FALSE])) + mean(colMeans(X[, 4:ncol(X), with = FALSE]))
    X_demean %*% t(X_demean)
})) / (nrow(X) - no - np + nr)
Sys.time() - start.time