Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
example of a test design and mixed effect analysis for measuring polywell performance and characterizing symmetry
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.