1aov <- function(formula, ...){ 2 UseMethod("aov") 3} 4aov.default <- stats::aov 5aov.design <- function (formula, ..., response = NULL, degree = NULL, FUN = mean, 6 use.center = FALSE) 7{ 8 if (!"design" %in% class(formula)) 9 stop("aov.design works on class design objects only") 10 di <- design.info(formula) 11 ## capture designs from conf.design 12 if (is.null(di)) 13 if (is.null(di)) stop("aov.design does not work for class design from package conf.design") 14 15 fo <- formula(formula, ..., response = response, degree = degree, 16 FUN = deparse(substitute(FUN)), use.center = use.center) 17 if (di$repeat.only | (length(grep("param", di$type)) > 0 & 18 length(grep("wide", di$type)) == 0) | (length(grep("center", 19 di$type)) > 0 & !use.center)) 20 aus <- aov(fo, data = model.frame(fo, data = NULL), ...) 21 else aus <- aov(fo, data = model.frame(fo, data = formula), 22 ...) 23 class(aus) <- c("aov.design", class(aus)) 24 25 if (length(grep("splitplot", di$type)) > 0) { 26 mm <- model.matrix(aus) 27 coefs <- coef(aus)[-1] 28 nms <- names(coefs) 29 hilf <- apply(mm[, 1 + (1:di$nfac.WP), drop = FALSE], 30 1, paste, collapse = "") 31 pchs <- rep("*", length(coefs)) 32 pchs[1:di$nfac.WP] <- "o" 33 for (j in setdiff(1:(length(coefs)), 1:di$nfac.WP)) { 34 if (!length(table(paste(hilf, mm[, nms[j]], sep = ""))) > 35 di$nWPs) 36 pchs[j] <- "o" 37 } 38 aus$WholePlotEffects <- names(coefs[pchs=="o"]) 39 } 40 aus 41} 42 43summary.aov.design <- function(object, ...){ 44 aus <- summary.aov(object) 45 class(aus) <- c("summary.aov.design", class(aus)) 46 fop <- formula(terms(object)) 47 attributes(fop) <- NULL 48 attr(aus, "formula") <- fop 49 attr(aus, "WholePlotEffects") <- object$WholePlotEffects 50 aus 51} 52 53print.summary.aov.design <- function(x, ...){ 54 cat("Number of observations used:", sum(x[[1]]$Df) + 1,"\n") 55 cat("Formula:\n") 56 print(attr(x, "formula")) 57 58 ## printing x with the summary.aov method 59 NextMethod(x, ...) 60 61 if (!is.null(attr(x, "WholePlotEffects"))){ 62 cat("WARNING: This is a split plot design, whole plot effects may have larger variance!\n") 63 cat(" p-values for whole plot effects may be misleadingly low!\n") 64 cat("The whole plot effects are:\n") 65 print(attr(x, "WholePlotEffects"), quote=FALSE) 66 }} 67 68print.aov.design <-function(x, ...){ 69 cat("Number of observations used:", nrow(model.frame(x)),"\n") 70 cat("Formula:\n") 71 fop <- formula(x) 72 attributes(fop) <- NULL 73 print(fop) 74 ## printing x with the aov method 75 NextMethod(x, ...) 76 77 if (!is.null(x$WholePlotEffects)){ 78 cat("WARNING: This is a split plot design, whole plot effects may have larger variance!\n") 79 cat(" p-values for whole plot effects may be misleadingly low!\n") 80 cat("The whole plot effects are:\n") 81 print(x$WholePlotEffects, quote=FALSE) 82 }} 83