library(preprocessCore)

###
### Creating a simulation model
###

col.effects <- c(10,11,10.5,12,9.5)
row.effects <- c(seq(-0.5,-0.1,by=0.1),seq(0.1,0.5,by=0.1))
 
y <- outer(row.effects, col.effects,"+")


###
### Checking that the estimates agree
###
### Note that this function fits
###   y = col + row + error
### with sum(row) = 0 
###

rcModelPLM(y)$Estimates
rcModelPLM(t(y))$Estimates

transposedEstimates <- rcModelPLM(t(y))$Estimates

transposedEstimates[1:nrow(y)] - mean(transposedEstimates[1:nrow(y)])
transposedEstimates[(nrow(y) + 1):(nrow(y) +ncol(y))] + mean(transposedEstimates[1:nrow(y)])


###
### Checking that residuals and Weights agree
###

w <- runif(50)
y <- y + rnorm(50)

rcModelPLM(y)$Estimates
rcModelPLM(t(y))$Estimates
transposedEstimates <- rcModelPLM(t(y))$Estimates
transposedEstimates[1:nrow(y)] - mean(transposedEstimates[1:nrow(y)])
transposedEstimates[(nrow(y) + 1):(nrow(y) +ncol(y))] + mean(transposedEstimates[1:nrow(y)])

###
### These should be at the level of machine error
###

max(abs(rcModelPLM(t(y))$Weights - t(rcModelPLM(y)$Weights)))
max(abs(rcModelPLM(t(y))$Residuals - t(rcModelPLM(y)$Residuals)))

###
### Unfortunately this is where things start breaking down
###

rcModelPLM(y)$StdErrors
rcModelPLM(t(y))$StdErrors


###
### Super-secret testing functions.
### Forget you ever saw them
### these are non-public
###
### Although this is not really my favorite
### way of computing the SE estimates, but
### it is the current default
###

residSE <- sqrt(sum(rcModelPLM(y)$Residuals^2*rcModelPLM(y)$Weights)/(50 - 14))

weights <- rcModelPLM(y)$Weights
XTWX <- .C("XTWX_R",as.integer(10),as.integer(5),as.double(weights),double(14^2),PACKAGE="preprocessCore")[[4]]
XTWXinv <- matrix(.C("XTWX_R_inv",as.integer(10),as.integer(5),as.double(XTWX))[[3]],14,14)


weights2 <- rcModelPLM(t(y))$Weights
XTWX2 <- .C("XTWX_R",as.integer(5),as.integer(10),as.double(weights2),double(14^2),PACKAGE="preprocessCore")[[4]]
XTWXinv2 <- matrix(.C("XTWX_R_inv",as.integer(5),as.integer(10),as.double(XTWX2))[[3]],14,14)


residSE*sqrt(diag(XTWXinv))
residSE*sqrt(diag(XTWXinv2))

### Compare with

rcModelPLM(y)$StdError
rcModelPLM(t(y))$StdError