example of a test design and mixed effect analysis for measuring polywell performance and characterizing symmetry

polywell_sym.R

# file: polywell_sym.R
#   example of a test design and mixed effect analysis for measuring
#   polywell performance and characterizing symmetry
# author: Josh Stults, www.variousconsequences.com 
# date: 28 July 2012
# see also:
#   http://prometheusfusionperfection.com/2012/07/27/symmetry-test/ 
#   http://cran.r-project.org/web/packages/AlgDesign/index.html
#   www.stat.wisc.edu/~bates/IMPS2008/lme4D.pdf

# clear the workspace: 
rm(list=ls())
# for doing design of experiments: 
library(AlgDesign)
# for fitting mixed effects models: 
library(lme4)

# 7 probe locations: six faces of the cube and the center 
probeloc = unique(rbind(gen.factorial(levels=c(3,1,1),nVars=3,varNames=c("x","y","z")),gen.factorial(levels=c(1,3,1),nVars=3,varNames=c("x","y","z")),gen.factorial(levels=c(1,1,3),nVars=3,varNames=c("x","y","z"))))

# three run design: 
dx.1 = optBlock(~x+y+z, withinData=probeloc, blocksizes=rep(3,3), wholeBlockData=data.frame(run=seq(1,3)))
# four run design:
dx.2 = optBlock(~x+y+z, withinData=probeloc, blocksizes=rep(3,4), wholeBlockData=data.frame(run=seq(1,4)))
# five run design: 
dx.3 = optBlock(~x+y+z, withinData=probeloc, blocksizes=rep(3,5), wholeBlockData=data.frame(run=seq(1,5)))

# fitting the mixed effects model to a dummy response:  
dummydat = transform(dx.3$design, dummy.response=rnorm(nrow(dx.3$design)))
lm.1 = lmer(dummy.response ~ x + y + z + (1|run), data=dummydat)