require("reshape2")
dat <- yates.missing
mat0 <- acast(dat[, c('trt','block','y')], trt~block,
id.var=c('trt','block'))
# Use lm to estimate missing values. The estimated missing values
# are the same as in Yates (1933)
m1 <- lm(y~trt+block, dat)
dat$pred <- predict(m1, new=dat[, c('trt','block')])
dat$filled <- ifelse(is.na(dat$y), dat$pred, dat$y)
mat1 <- acast(dat[, c('trt','block','pred')], trt~block,
id.var=c('trt','block'))
# Another method to estimate missing values via PCA
require("pcaMethods")
m2 <- pca(mat0, method="nipals", center=FALSE, nPcs=3)
mat2 <- m2@scores
# Compare
ord <- c("O","N","K","P","NK","NP","KP","NKP")
print(mat0[ord,], na.print=".")
round(mat1[ord,] ,2)
round(mat2[ord,] ,2)
# SVD with 3 components recovers original data better
sum((mat0-mat1)^2, na.rm=TRUE)
sum((mat0-mat2)^2, na.rm=TRUE) # Smaller SS => better fitRun the code above in your browser using DataLab