#'# Script to calculate power by simulation for the Farewell and Herzberg example #' #'## Initialize library(knitr) #knitr::spin("FHpower.r") library(dae) packageVersion("dae") library(asreml) packageVersion("asreml") library(asremlPlus) packageVersion("asremlPlus") library(parallel) library(foreach) library(doParallel) options(width = 95) source("globalPower.r") asreml.options(ai.sing = TRUE) #nsimul <- 5 #'## Load the design load("plaid.dat.rda") #'## Add factors for heterogeneous variances plaid.dat <- within(plaid.dat, { MotionsExpress <- fac.combine(list(Motions, Expressiveness), combine.levels = TRUE) WMotionsExpress <- fac.nested(MotionsExpress) }) plaid.dat <- with(plaid.dat, plaid.dat[order(MotionsExpress, WMotionsExpress), ]) sim.dat <- plaid.dat[, -match("Y", names(plaid.dat))] n <- nrow(plaid.dat) #'## Generate the simulation data sets FH.dat <- makeSimData(params = params, design = sim.dat, TME = TME, nsimul = nsimul) #'## Initialize parallel processing cl <- makeCluster(ncores) registerDoParallel(cl) cat("\n\n#### Number of simulated data sets is ",nsimul,"\n") #'## Analyze the simulated data for model d FH.d.stats <- foreach(k = 1:nsimul, .errorhandling = "pass", .inorder=TRUE, .packages = c("asreml", "asremlPlus", "dae")) %dopar% { analPower(k, dat = FH.dat$data) } FH.d.stats <- do.call(rbind, FH.d.stats) head(FH.d.stats, n = 25) summary(FH.d.stats) #'## Analyze the simulated data for model c FH.c.stats <- foreach(k = 1:nsimul, .errorhandling = "pass", .inorder=TRUE, .packages = c("asreml", "asremlPlus", "dae")) %dopar% { analPower(k, dat = FH.dat$data, model = "c") } FH.c.stats <- do.call(rbind, FH.c.stats) head(FH.c.stats, n = 25) summary(FH.c.stats) #'## Hand back clusters stopCluster(cl) #'## Calculate power statistics getStats <- function(stats) { converged <- table(stats$ModelFit)["F"] power <- sum(stats$pTME <= 0.05, na.rm = TRUE)/as.numeric(table(stats$ModelFit)["F"]) return(cbind(converged, power)) } #'## Power results cat("\n\n#### Power statistics for model c)\n\n") getStats(FH.c.stats) cat("\n\n#### Power statistics for model d)\n\n") getStats(FH.d.stats) #'## Save the workspace save.image("FHpower.RData")