R version 4.2.1 (2022-06-23) -- "Funny-Looking Kid" Copyright (C) 2022 The R Foundation for Statistical Computing Platform: aarch64-apple-darwin20 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > pkgname <- "afex" > source(file.path(R.home("share"), "R", "examples-header.R")) > options(warn = 1) > library('afex') Loading required package: lme4 Loading required package: Matrix ************ Welcome to afex. For support visit: http://afex.singmann.science/ - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4() - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB' - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default. - Get and set global package options with: afex_options() - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts() - For example analyses see: browseVignettes("afex") ************ Attaching package: ‘afex’ The following object is masked from ‘package:lme4’: lmer > > base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') > base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv') > cleanEx() > nameEx("afex_options") > ### * afex_options > > flush(stderr()); flush(stdout()) > > ### Name: afex_options > ### Title: Set/get global afex options > ### Aliases: afex_options > > ### ** Examples > > afex_options() # see all options $check_contrasts [1] TRUE $correction_aov [1] "GG" $emmeans_model [1] "multivariate" $es_aov [1] "ges" $factorize [1] TRUE $include_aov [1] FALSE $lmer_function [1] "lmerTest" $method_mixed [1] "S" $return_aov [1] "afex_aov" $round_ps function (x) { substr(as.character(ifelse(x < 0.001, " <.001", ifelse(round(x, 3) == 1, " >.999", formatC(x, digits = 3, format = "f")))), 2, 7) } $set_data_arg [1] FALSE $sig_symbols [1] " +" " *" " **" " ***" $type [1] 3 > > afex_options("return_aov") #get single option [1] "afex_aov" > > aop <- afex_options() # save current options > > ## Not run: > ##D # change options > ##D afex_options(return_aov = "nice") > ##D afex_options("return_aov") #get single option > ##D afex_options(return_aov = "nice", method_mixed = "LRT") > ##D afex_options("method_mixed") #get single option > ##D # do something > ## End(Not run) > afex_options(aop) # reset options > > > > > cleanEx() > nameEx("afex_plot") > ### * afex_plot > > flush(stderr()); flush(stdout()) > > ### Name: afex_plot > ### Title: m-way Plot with Error Bars and Raw Data > ### Aliases: afex_plot afex_plot.afex_aov afex_plot.mixed afex_plot.merMod > ### afex_plot.default interaction_plot oneway_plot > > ### ** Examples > > # note: use library("ggplot") to avoid "ggplot2::" in the following > > ################################################################## > ## 2-factor Within-Subject Design ## > ################################################################## > > data(md_12.1) > aw <- aov_ez("id", "rt", md_12.1, within = c("angle", "noise")) > > ##--------------------------------------------------------------- > ## Basic Interaction Plots - > ##--------------------------------------------------------------- > > ## all examples require emmeans and ggplot2: > if (requireNamespace("emmeans") && requireNamespace("ggplot2")) { + + afex_plot(aw, x = "angle", trace = "noise") + # or: afex_plot(aw, x = ~angle, trace = ~noise) + + afex_plot(aw, x = "noise", trace = "angle") + + ### For within-subject designs, using within-subject CIs is better: + afex_plot(aw, x = "angle", trace = "noise", error = "within") + (p1 <- afex_plot(aw, x = "noise", trace = "angle", error = "within")) + + ## use different themes for nicer graphs: + p1 + ggplot2::theme_bw() + } Warning: Panel(s) show within-subjects factors, but not within-subjects error bars. For within-subjects error bars use: error = "within" Warning: Panel(s) show within-subjects factors, but not within-subjects error bars. For within-subjects error bars use: error = "within" > ## Not run: > ##D p1 + ggplot2::theme_light() > ##D p1 + ggplot2::theme_minimal() > ##D p1 + jtools::theme_apa() > ##D p1 + ggpubr::theme_pubr() > ##D > ##D ### set theme globally for R session: > ##D ggplot2::theme_set(ggplot2::theme_bw()) > ##D > ##D ### There are several ways to deal with overlapping points in the background besides alpha > ##D # 1. using the default data geom and ggplot2::position_jitterdodge > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", dodge = 0.3, > ##D data_arg = list( > ##D position = > ##D ggplot2::position_jitterdodge( > ##D jitter.width = 0, > ##D jitter.height = 5, > ##D dodge.width = 0.3 ## needs to be same as dodge > ##D ), > ##D color = "darkgrey")) > ##D > ##D # 2. using ggbeeswarm::geom_beeswarm > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", dodge = 0.5, > ##D data_geom = ggbeeswarm::geom_beeswarm, > ##D data_arg = list( > ##D dodge.width = 0.5, ## needs to be same as dodge > ##D cex = 0.8, > ##D color = "darkgrey")) > ##D > ##D # 3. do not display points, but use a violinplot: ggplot2::geom_violin > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", > ##D data_geom = ggplot2::geom_violin, > ##D data_arg = list(width = 0.5)) > ##D > ##D # 4. violinplots with color: ggplot2::geom_violin > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", > ##D mapping = c("linetype", "shape", "fill"), > ##D data_geom = ggplot2::geom_violin, > ##D data_arg = list(width = 0.5)) > ##D > ##D # 5. do not display points, but use a boxplot: ggplot2::geom_boxplot > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", > ##D data_geom = ggplot2::geom_boxplot, > ##D data_arg = list(width = 0.3)) > ##D > ##D # 6. combine points with boxplot: ggpol::geom_boxjitter > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", > ##D data_geom = ggpol::geom_boxjitter, > ##D data_arg = list(width = 0.3)) > ##D ## hides error bars! > ##D > ##D # 7. nicer variant of ggpol::geom_boxjitter > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", > ##D mapping = c("shape", "fill"), > ##D data_geom = ggpol::geom_boxjitter, > ##D data_arg = list( > ##D width = 0.3, > ##D jitter.width = 0, > ##D jitter.height = 10, > ##D outlier.intersect = TRUE), > ##D point_arg = list(size = 2.5), > ##D error_arg = list(size = 1.5, width = 0)) > ##D > ##D # 8. nicer variant of ggpol::geom_boxjitter without lines > ##D afex_plot(aw, x = "noise", trace = "angle", error = "within", dodge = 0.7, > ##D mapping = c("shape", "fill"), > ##D data_geom = ggpol::geom_boxjitter, > ##D data_arg = list( > ##D width = 0.5, > ##D jitter.width = 0, > ##D jitter.height = 10, > ##D outlier.intersect = TRUE), > ##D point_arg = list(size = 2.5), > ##D line_arg = list(linetype = 0), > ##D error_arg = list(size = 1.5, width = 0)) > ## End(Not run) > > > ##--------------------------------------------------------------- > ## One-Way Plots - > ##--------------------------------------------------------------- > > ## Not run: > ##D afex_plot(aw, x = "angle", error = "within") ## default > ##D > ##D ## with color we need larger points > ##D afex_plot(aw, x = "angle", mapping = "color", error = "within", > ##D point_arg = list(size = 2.5), > ##D error_arg = list(size = 1.5, width = 0.05)) > ##D > ##D afex_plot(aw, x = "angle", error = "within", data_geom = ggpol::geom_boxjitter) > ##D > ##D ## nicer > ##D afex_plot(aw, x = "angle", error = "within", data_geom = ggpol::geom_boxjitter, > ##D mapping = "fill", data_alpha = 0.7, > ##D data_arg = list( > ##D width = 0.6, > ##D jitter.width = 0.07, > ##D jitter.height = 10, > ##D outlier.intersect = TRUE > ##D ), > ##D point_arg = list(size = 2.5), > ##D error_arg = list(size = 1.5, width = 0.05)) > ##D > ##D ## we can add a line connecting the means using geom_point(aes(group = 1)): > ##D afex_plot(aw, x = "angle", error = "within") + > ##D ggplot2::geom_line(ggplot2::aes(group = 1)) > ##D > ##D ## One-way plots also supports panels: > ##D afex_plot(aw, x = "angle", panel = "noise", error = "within") > ##D > ##D ## And panels with lines: > ##D afex_plot(aw, x = "angle", panel = "noise", error = "within") + > ##D ggplot2::geom_line(ggplot2::aes(group = 1)) > ##D > ##D > ##D ## For more complicated plots it is easier to attach ggplot2: > ##D library("ggplot2") > ##D > ##D ## We can hide geoms by plotting them in transparent color and add them > ##D ## afterward to use a mapping not directly supported. > ##D ## For example, the next plot adds a line to a one-way plot with panels, but > ##D ## with all geoms in the foreground having a color conditional on the panel. > ##D > ##D afex_plot(aw, x = "angle", panel = "noise", error = "within", > ##D point_arg = list(color = "transparent"), > ##D error_arg = list(color = "transparent")) + > ##D geom_point(aes(color = panel)) + > ##D geom_linerange(aes(color = panel, ymin = lower, ymax = upper)) + > ##D geom_line(aes(group = 1, color = panel)) + > ##D guides(color = guide_legend(title = "NOISE")) > ##D ## Note that we need to use guides explicitly, otherwise the legend title would > ##D ## be "panel". legend_title does not work in this case. > ##D > ##D ##--------------------------------------------------------------- > ##D ## Other Basic Options - > ##D ##--------------------------------------------------------------- > ##D > ##D ## relabel factor levels via factor_levels (with message) > ##D afex_plot(aw, x = "noise", trace = "angle", > ##D factor_levels = list(angle = c("0°", "4°", "8°"), > ##D noise = c("Absent", "Present"))) > ##D > ##D ## factor_levels allows named vectors which enable reordering the factor levels > ##D ### and renaming subsets of levels: > ##D afex_plot(aw, x = "noise", trace = "angle", > ##D factor_levels = list( > ##D angle = c(X8 = "8°", X4 = "4°", X0 = "0°"), > ##D noise = c(present = "Present") > ##D ) > ##D ) > ##D > ##D > ##D ## Change title of legend > ##D afex_plot(aw, x = "noise", trace = "angle", > ##D legend_title = "Noise Condition") > ##D > ##D ## for plots with few factor levels, smaller dodge might be better: > ##D afex_plot(aw, x = "angle", trace = "noise", dodge = 0.25) > ##D > ##D ################################################################# > ##D ## 4-factor Mixed Design ## > ##D ################################################################# > ##D > ##D data(obk.long, package = "afex") > ##D a1 <- aov_car(value ~ treatment * gender + Error(id/(phase*hour)), > ##D data = obk.long, observed = "gender") > ##D > ##D ## too difficult to see anything > ##D afex_plot(a1, ~phase*hour, ~treatment) + > ##D ggplot2::theme_light() > ##D > ##D ## better > ##D afex_plot(a1, ~hour, ~treatment, ~phase) + > ##D ggplot2::theme_light() > ##D > ##D ## even better > ##D afex_plot(a1, ~hour, ~treatment, ~phase, > ##D dodge = 0.65, > ##D data_arg = list( > ##D position = > ##D ggplot2::position_jitterdodge( > ##D jitter.width = 0, > ##D jitter.height = 0.2, > ##D dodge.width = 0.65 ## needs to be same as dodge > ##D ), > ##D color = "darkgrey")) + > ##D ggplot2::theme_classic() > ##D > ##D # with color instead of linetype to separate trace factor > ##D afex_plot(a1, ~hour, ~treatment, ~phase, > ##D mapping = c("shape", "color"), > ##D dodge = 0.65, > ##D data_arg = list( > ##D position = > ##D ggplot2::position_jitterdodge( > ##D jitter.width = 0, > ##D jitter.height = 0.2, > ##D dodge.width = 0.65 ## needs to be same as dodge > ##D ))) + > ##D ggplot2::theme_light() > ##D > ##D # only color to separate trace factor > ##D afex_plot(a1, ~hour, ~treatment, ~phase, > ##D mapping = "color", > ##D dodge = 0.65, > ##D data_arg = list( > ##D position = > ##D ggplot2::position_jitterdodge( > ##D jitter.width = 0, > ##D jitter.height = 0.2, > ##D dodge.width = 0.65 ## needs to be same as dodge > ##D ))) + > ##D ggplot2::theme_classic() > ##D > ##D > ##D ## plot involving all 4 factors: > ##D afex_plot(a1, ~hour, ~treatment, ~gender+phase, > ##D dodge = 0.65, > ##D data_arg = list( > ##D position = > ##D ggplot2::position_jitterdodge( > ##D jitter.width = 0, > ##D jitter.height = 0.2, > ##D dodge.width = 0.65 ## needs to be same as dodge > ##D ), > ##D color = "darkgrey")) + > ##D ggplot2::theme_bw() > ##D > ##D > ##D ##--------------------------------------------------------------- > ##D ## Different Standard Errors Available - > ##D ##--------------------------------------------------------------- > ##D > ##D ## purely within-design > ##D cbind( > ##D afex_plot(a1, ~phase, ~hour, > ##D error = "model", return = "data")$means[,c("phase", "hour", "y", "SE")], > ##D multivariate = afex_plot(a1, ~phase, ~hour, > ##D error = "model", return = "data")$means$error, > ##D mean = afex_plot(a1, ~phase, ~hour, > ##D error = "mean", return = "data")$means$error, > ##D within = afex_plot(a1, ~phase, ~hour, > ##D error = "within", return = "data")$means$error, > ##D between = afex_plot(a1, ~phase, ~hour, > ##D error = "between", return = "data")$means$error) > ##D ## mixed design > ##D cbind( > ##D afex_plot(a1, ~phase, ~treatment, > ##D error = "model", return = "data")$means[,c("phase", "treatment", "y", "SE")], > ##D multivariate = afex_plot(a1, ~phase, ~treatment, > ##D error = "model", return = "data")$means$error, > ##D mean = afex_plot(a1, ~phase, ~treatment, > ##D error = "mean", return = "data")$means$error, > ##D within = afex_plot(a1, ~phase, ~treatment, > ##D error = "within", return = "data")$means$error, > ##D between = afex_plot(a1, ~phase, ~treatment, > ##D error = "between", return = "data")$means$error) > ## End(Not run) > > ################################################################## > ## Mixed Models ## > ################################################################## > if (requireNamespace("MEMSS") && + requireNamespace("emmeans") && + requireNamespace("ggplot2")) { + + data("Machines", package = "MEMSS") + m1 <- mixed(score ~ Machine + (Machine|Worker), data=Machines) + pairs(emmeans::emmeans(m1, "Machine")) + # contrast estimate SE df t.ratio p.value + # A - B -7.966667 2.420850 5 -3.291 0.0481 + # A - C -13.916667 1.540100 5 -9.036 0.0007 + # B - C -5.950000 2.446475 5 -2.432 0.1253 + + ## Default (i.e., model-based) error bars suggest no difference between Machines. + ## This contrasts with pairwise comparisons above. + afex_plot(m1, "Machine") + + ## Impression from within-subject error bars is more in line with pattern of differences. + afex_plot(m1, "Machine", error = "within") + } Loading required namespace: MEMSS Contrasts set to contr.sum for the following variables: Machine, Worker Aggregating data over: Worker Aggregating data over: Worker > > ## Not run: > ##D data("fhch2010") # load > ##D fhch <- droplevels(fhch2010[ fhch2010$correct,]) # remove errors > ##D ### following model should take less than a minute to fit: > ##D mrt <- mixed(log_rt ~ task*stimulus*frequency + (stimulus*frequency||id)+ > ##D (task||item), fhch, method = "S", expand_re = TRUE) > ##D > ##D ## way too many points in background: > ##D afex_plot(mrt, "stimulus", "frequency", "task") > ##D > ##D ## better to restrict plot of data to one random-effects grouping variable > ##D afex_plot(mrt, "stimulus", "frequency", "task", id = "id") > ##D ## when plotting data from a single random effect, different error bars are possible: > ##D afex_plot(mrt, "stimulus", "frequency", "task", id = "id", error = "within") > ##D afex_plot(mrt, "stimulus", "frequency", "task", id = "id", error = "mean") > ##D > ##D ## compare visual impression with: > ##D pairs(emmeans::emmeans(mrt, c("stimulus", "frequency"), by = "task")) > ##D > ##D ## same logic also possible for other random-effects grouping factor > ##D afex_plot(mrt, "stimulus", "frequency", "task", id = "item") > ##D ## within-item error bars are misleading here. task is sole within-items factor. > ##D afex_plot(mrt, "stimulus", "frequency", "task", id = "item", error = "within") > ##D ## CIs based on stanard error of mean look small, but not unreasonable given results. > ##D afex_plot(mrt, "stimulus", "frequency", "task", id = "item", error = "mean") > ##D > ##D ### compare distribution of individual data for different random effects: > ##D ## requires package cowplot > ##D p_id <- afex_plot(mrt, "stimulus", "frequency", "task", id = "id", > ##D error = "within", dodge = 0.7, > ##D data_geom = ggplot2::geom_violin, > ##D mapping = c("shape", "fill"), > ##D data_arg = list(width = 0.7)) + > ##D ggplot2::scale_shape_manual(values = c(4, 17)) + > ##D ggplot2::labs(title = "ID") > ##D > ##D p_item <- afex_plot(mrt, "stimulus", "frequency", "task", id = "item", > ##D error = "within", dodge = 0.7, > ##D data_geom = ggplot2::geom_violin, > ##D mapping = c("shape", "fill"), > ##D data_arg = list(width = 0.7)) + > ##D ggplot2::scale_shape_manual(values = c(4, 17)) + > ##D ggplot2::labs(title = "Item") > ##D > ##D ### see: https://cran.r-project.org/package=cowplot/vignettes/shared_legends.html > ##D p_comb <- cowplot::plot_grid( > ##D p_id + ggplot2::theme_light() + ggplot2::theme(legend.position="none"), > ##D p_item + ggplot2::theme_light() + ggplot2::theme(legend.position="none") > ##D ) > ##D legend <- cowplot::get_legend(p_id + ggplot2::theme(legend.position="bottom")) > ##D cowplot::plot_grid(p_comb, legend, > ##D ncol = 1, > ##D rel_heights = c(1, 0.1)) > ##D > ##D ##---------------------------------------------------------------- > ##D ## Support for lme4::lmer - > ##D ##---------------------------------------------------------------- > ##D > ##D Oats <- nlme::Oats > ##D ## afex_plot does currently not support implicit nesting: (1|Block/Variety) > ##D ## Instead, we need to create the factor explicitly > ##D Oats$VarBlock <- Oats$Variety:Oats$Block > ##D Oats.lmer <- lmer(yield ~ Variety * factor(nitro) + (1|VarBlock) + (1|Block), > ##D data = Oats) > ##D afex_plot(Oats.lmer, "nitro", "Variety") > ##D afex_plot(Oats.lmer, "nitro", panel = "Variety") > ##D > ##D ################################################################## > ##D ## Default Method works for Models Supported by emmeans ## > ##D ################################################################## > ##D > ##D ## lm > ##D warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks) > ##D afex_plot(warp.lm, "tension") > ##D afex_plot(warp.lm, "tension", "wool") > ##D > ##D ## poisson glm > ##D ins <- data.frame( > ##D n = c(500, 1200, 100, 400, 500, 300), > ##D size = factor(rep(1:3,2), labels = c("S","M","L")), > ##D age = factor(rep(1:2, each = 3)), > ##D claims = c(42, 37, 1, 101, 73, 14)) > ##D ins.glm <- glm(claims ~ size + age + offset(log(n)), > ##D data = ins, family = "poisson") > ##D afex_plot(ins.glm, "size", "age") > ##D > ##D ## binomial glm adapted from ?predict.glm > ##D ldose <- factor(rep(0:5, 2)) > ##D numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16) > ##D sex <- factor(rep(c("M", "F"), c(6, 6))) > ##D SF <- numdead/20 ## dv should be a vector, no matrix > ##D budworm.lg <- glm(SF ~ sex*ldose, family = binomial, > ##D weights = rep(20, length(numdead))) > ##D afex_plot(budworm.lg, "ldose") > ##D afex_plot(budworm.lg, "ldose", "sex") ## data point is hidden behind mean! > ##D afex_plot(budworm.lg, "ldose", "sex", > ##D data_arg = list(size = 4, color = "red")) > ##D > ##D ## nlme mixed model > ##D data(Oats, package = "nlme") > ##D Oats$nitro <- factor(Oats$nitro) > ##D oats.1 <- nlme::lme(yield ~ nitro * Variety, > ##D random = ~ 1 | Block / Variety, > ##D data = Oats) > ##D afex_plot(oats.1, "nitro", "Variety", data = Oats) > ##D afex_plot(oats.1, "nitro", "Variety", data = Oats, id = "Block") > ##D afex_plot(oats.1, "nitro", data = Oats) > ##D afex_plot(oats.1, "nitro", data = Oats, id = c("Block", "Variety")) > ##D afex_plot(oats.1, "nitro", data = Oats, id = "Block") > ##D > ## End(Not run) > > > > cleanEx() > nameEx("all_fit") > ### * all_fit > > flush(stderr()); flush(stdout()) > > ### Name: all_fit > ### Title: Refit 'lmer' model using multiple optimizers > ### Aliases: all_fit nmkbw > > ### ** Examples > > > ## Not run: > ##D > ##D # basic usage > ##D require(optimx) > ##D gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), > ##D data = cbpp, family = binomial) > ##D gm_all <- all_fit(gm1) > ##D t(sapply(gm_all,fixef)) ## extract fixed effects > ##D sapply(gm_all,logLik) ## log-likelihoods > ##D sapply(gm_all,getME,"theta") ## theta parameters > ##D !sapply(gm_all,inherits,"try-error") ## was fit OK? > ##D > ##D > ##D ## for GLMMs: > ##D require("mlmRev") # for data > ##D gm1 <- mixed(use ~ age*urban + (1 | district), family = binomial, > ##D data = Contraception, method = "LRT") > ##D gm_all <- all_fit(gm1$full_model) > ##D sapply(gm_all,logLik) > ##D > ##D ## use allFit in combination with expand.re = TRUE > ##D data("sk2011.2") # see example("mixed") > ##D sk2_aff <- droplevels(sk2011.2[sk2011.2$what == "affirmation",]) > ##D sk_m2 <- mixed(response ~ instruction*inference*type+(inference*type||id), sk2_aff, > ##D expand_re = TRUE) > ##D sk_m2 > ##D sk_m2_allFit <- all_fit(sk_m2$full_model) > ##D sk_m2_allFit # all fits fail > ##D > ##D sk_m2_allFit <- all_fit(sk_m2$full_model, data = sk_m2$data) # works > ##D t(sapply(sk_m2_allFit,fixef)) > ##D sapply(sk_m2_allFit,logLik) > ##D > ## End(Not run) > > > > cleanEx() > nameEx("aov_car") > ### * aov_car > > flush(stderr()); flush(stdout()) > > ### Encoding: UTF-8 > > ### Name: aov_car > ### Title: Convenient ANOVA estimation for factorial designs > ### Aliases: aov_car aov_4 aov_ez > > ### ** Examples > > > ########################## > ## 1: Specifying ANOVAs ## > ########################## > > # Example using a purely within-subjects design > # (Maxwell & Delaney, 2004, Chapter 12, Table 12.5, p. 578): > data(md_12.1) > aov_ez("id", "rt", md_12.1, within = c("angle", "noise"), + anova_table=list(correction = "none", es = "none")) Anova Table (Type 3 tests) Response: rt Effect df MSE F p.value 1 angle 2, 18 3560.00 40.72 *** <.001 2 noise 1, 9 8460.00 33.77 *** <.001 3 angle:noise 2, 18 1160.00 45.31 *** <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > > # Default output > aov_ez("id", "rt", md_12.1, within = c("angle", "noise")) Anova Table (Type 3 tests) Response: rt Effect df MSE F ges p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .390 <.001 2 noise 1, 9 8460.00 33.77 *** .387 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .188 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > > # examples using obk.long (see ?obk.long), a long version of the OBrienKaiser > # dataset (car package). Data is a split-plot or mixed design: contains both > # within- and between-subjects factors. > data(obk.long, package = "afex") > > # estimate mixed ANOVA on the full design: > aov_car(value ~ treatment * gender + Error(id/(phase*hour)), + data = obk.long, observed = "gender") Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 22.81 3.94 + .198 .055 2 gender 1, 10 22.81 3.66 + .115 .085 3 treatment:gender 2, 10 22.81 2.86 .179 .104 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 8 hour 1.84, 18.41 3.39 16.69 *** .125 <.001 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_4(value ~ treatment * gender + (phase*hour|id), + data = obk.long, observed = "gender") Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 22.81 3.94 + .198 .055 2 gender 1, 10 22.81 3.66 + .115 .085 3 treatment:gender 2, 10 22.81 2.86 .179 .104 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 8 hour 1.84, 18.41 3.39 16.69 *** .125 <.001 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_ez("id", "value", obk.long, between = c("treatment", "gender"), + within = c("phase", "hour"), observed = "gender") Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 22.81 3.94 + .198 .055 2 gender 1, 10 22.81 3.66 + .115 .085 3 treatment:gender 2, 10 22.81 2.86 .179 .104 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 8 hour 1.84, 18.41 3.39 16.69 *** .125 <.001 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > # the three calls return the same ANOVA table: > # Anova Table (Type 3 tests) > # > # Response: value > # Effect df MSE F ges p.value > # 1 treatment 2, 10 22.81 3.94 + .198 .055 > # 2 gender 1, 10 22.81 3.66 + .115 .085 > # 3 treatment:gender 2, 10 22.81 2.86 .179 .104 > # 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001 > # 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013 > # 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 > # 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 > # 8 hour 1.84, 18.41 3.39 16.69 *** .125 <.001 > # 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979 > # 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628 > # 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641 > # 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335 > # 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930 > # 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449 > # 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646 > # --- > # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > # > # Sphericity correction method: GG > > # "numeric" variables are per default converted to factors (as long as factorize > # = TRUE): > obk.long$hour2 <- as.numeric(as.character(obk.long$hour)) > > # gives same results as calls before > aov_car(value ~ treatment * gender + Error(id/phase*hour2), + data = obk.long, observed = c("gender")) Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 22.81 3.94 + .198 .055 2 gender 1, 10 22.81 3.66 + .115 .085 3 treatment:gender 2, 10 22.81 2.86 .179 .104 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 8 hour2 1.84, 18.41 3.39 16.69 *** .125 <.001 9 treatment:hour2 3.68, 18.41 3.39 0.09 .002 .979 10 gender:hour2 1.84, 18.41 3.39 0.45 .004 .628 11 treatment:gender:hour2 3.68, 18.41 3.39 0.62 .011 .641 12 phase:hour2 3.60, 35.96 2.67 1.18 .015 .335 13 treatment:phase:hour2 7.19, 35.96 2.67 0.35 .009 .930 14 gender:phase:hour2 3.60, 35.96 2.67 0.93 .012 .449 15 treatment:gender:phase:hour2 7.19, 35.96 2.67 0.74 .019 .646 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > > # ANCOVA: adding a covariate (necessary to set factorize = FALSE) > aov_car(value ~ treatment * gender + age + Error(id/(phase*hour)), + data = obk.long, observed = c("gender", "age"), factorize = FALSE) Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 9 23.96 3.58 + .202 .072 2 gender 1, 9 23.96 3.95 + .140 .078 3 age 1, 9 23.96 0.52 .018 .490 4 treatment:gender 2, 9 23.96 1.28 .091 .323 5 phase 1.70, 15.28 3.91 20.28 *** .166 <.001 6 treatment:phase 3.39, 15.28 3.91 6.07 ** .106 .005 7 gender:phase 1.70, 15.28 3.91 0.25 .002 .749 8 age:phase 1.70, 15.28 3.91 3.10 + .030 .081 9 treatment:gender:phase 3.39, 15.28 3.91 1.60 .031 .228 10 hour 2.14, 19.23 2.48 20.52 *** .138 <.001 11 treatment:hour 4.27, 19.23 2.48 0.71 .011 .601 12 gender:hour 2.14, 19.23 2.48 0.71 .006 .514 13 age:hour 2.14, 19.23 2.48 2.82 + .022 .082 14 treatment:gender:hour 4.27, 19.23 2.48 0.59 .009 .684 15 phase:hour 3.48, 31.36 2.83 0.99 .014 .419 16 treatment:phase:hour 6.97, 31.36 2.83 0.33 .010 .932 17 gender:phase:hour 3.48, 31.36 2.83 0.90 .013 .465 18 age:phase:hour 3.48, 31.36 2.83 0.77 .011 .540 19 treatment:gender:phase:hour 6.97, 31.36 2.83 0.65 .019 .710 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_4(value ~ treatment * gender + age + (phase*hour|id), + data = obk.long, observed = c("gender", "age"), factorize = FALSE) Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 9 23.96 3.58 + .202 .072 2 gender 1, 9 23.96 3.95 + .140 .078 3 age 1, 9 23.96 0.52 .018 .490 4 treatment:gender 2, 9 23.96 1.28 .091 .323 5 phase 1.70, 15.28 3.91 20.28 *** .166 <.001 6 treatment:phase 3.39, 15.28 3.91 6.07 ** .106 .005 7 gender:phase 1.70, 15.28 3.91 0.25 .002 .749 8 age:phase 1.70, 15.28 3.91 3.10 + .030 .081 9 treatment:gender:phase 3.39, 15.28 3.91 1.60 .031 .228 10 hour 2.14, 19.23 2.48 20.52 *** .138 <.001 11 treatment:hour 4.27, 19.23 2.48 0.71 .011 .601 12 gender:hour 2.14, 19.23 2.48 0.71 .006 .514 13 age:hour 2.14, 19.23 2.48 2.82 + .022 .082 14 treatment:gender:hour 4.27, 19.23 2.48 0.59 .009 .684 15 phase:hour 3.48, 31.36 2.83 0.99 .014 .419 16 treatment:phase:hour 6.97, 31.36 2.83 0.33 .010 .932 17 gender:phase:hour 3.48, 31.36 2.83 0.90 .013 .465 18 age:phase:hour 3.48, 31.36 2.83 0.77 .011 .540 19 treatment:gender:phase:hour 6.97, 31.36 2.83 0.65 .019 .710 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_ez("id", "value", obk.long, between = c("treatment", "gender"), + within = c("phase", "hour"), covariate = "age", + observed = c("gender", "age"), factorize = FALSE) Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 9 23.96 3.58 + .202 .072 2 gender 1, 9 23.96 3.95 + .140 .078 3 age 1, 9 23.96 0.52 .018 .490 4 treatment:gender 2, 9 23.96 1.28 .091 .323 5 phase 1.70, 15.28 3.91 20.28 *** .166 <.001 6 treatment:phase 3.39, 15.28 3.91 6.07 ** .106 .005 7 gender:phase 1.70, 15.28 3.91 0.25 .002 .749 8 age:phase 1.70, 15.28 3.91 3.10 + .030 .081 9 treatment:gender:phase 3.39, 15.28 3.91 1.60 .031 .228 10 hour 2.14, 19.23 2.48 20.52 *** .138 <.001 11 treatment:hour 4.27, 19.23 2.48 0.71 .011 .601 12 gender:hour 2.14, 19.23 2.48 0.71 .006 .514 13 age:hour 2.14, 19.23 2.48 2.82 + .022 .082 14 treatment:gender:hour 4.27, 19.23 2.48 0.59 .009 .684 15 phase:hour 3.48, 31.36 2.83 0.99 .014 .419 16 treatment:phase:hour 6.97, 31.36 2.83 0.33 .010 .932 17 gender:phase:hour 3.48, 31.36 2.83 0.90 .013 .465 18 age:phase:hour 3.48, 31.36 2.83 0.77 .011 .540 19 treatment:gender:phase:hour 6.97, 31.36 2.83 0.65 .019 .710 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > > # aggregating over one within-subjects factor (phase), with warning: > aov_car(value ~ treatment * gender + Error(id/hour), data = obk.long, + observed = "gender") Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 7.60 3.94 + .259 .055 2 gender 1, 10 7.60 3.66 + .162 .085 3 treatment:gender 2, 10 7.60 2.86 .253 .104 4 hour 1.84, 18.41 1.13 16.69 *** .168 <.001 5 treatment:hour 3.68, 18.41 1.13 0.09 .002 .979 6 gender:hour 1.84, 18.41 1.13 0.45 .005 .628 7 treatment:gender:hour 3.68, 18.41 1.13 0.62 .015 .641 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_ez("id", "value", obk.long, c("treatment", "gender"), "hour", + observed = "gender") Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 7.60 3.94 + .259 .055 2 gender 1, 10 7.60 3.66 + .162 .085 3 treatment:gender 2, 10 7.60 2.86 .253 .104 4 hour 1.84, 18.41 1.13 16.69 *** .168 <.001 5 treatment:hour 3.68, 18.41 1.13 0.09 .002 .979 6 gender:hour 1.84, 18.41 1.13 0.45 .005 .628 7 treatment:gender:hour 3.68, 18.41 1.13 0.62 .015 .641 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > # aggregating over both within-subjects factors (again with warning), > # only between-subjects factors: > aov_car(value ~ treatment * gender + Error(id), data = obk.long, + observed = c("gender")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 1.52 3.94 + .289 .055 2 gender 1, 10 1.52 3.66 + .189 .085 3 treatment:gender 2, 10 1.52 2.86 .295 .104 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > aov_4(value ~ treatment * gender + (1|id), data = obk.long, + observed = c("gender")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 1.52 3.94 + .289 .055 2 gender 1, 10 1.52 3.66 + .189 .085 3 treatment:gender 2, 10 1.52 2.86 .295 .104 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > aov_ez("id", "value", obk.long, between = c("treatment", "gender"), + observed = "gender") Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 1.52 3.94 + .289 .055 2 gender 1, 10 1.52 3.66 + .189 .085 3 treatment:gender 2, 10 1.52 2.86 .295 .104 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > > # only within-subject factors (ignoring between-subjects factors) > aov_car(value ~ Error(id/(phase*hour)), data = obk.long) Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 phase 1.54, 23.16 7.30 14.85 *** .147 <.001 2 hour 1.99, 29.91 2.46 21.63 *** .099 <.001 3 phase:hour 4.10, 61.56 2.00 1.35 .011 .260 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > aov_4(value ~ (phase*hour|id), data = obk.long) Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 phase 1.54, 23.16 7.30 14.85 *** .147 <.001 2 hour 1.99, 29.91 2.46 21.63 *** .099 <.001 3 phase:hour 4.10, 61.56 2.00 1.35 .011 .260 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > aov_ez("id", "value", obk.long, within = c("phase", "hour")) Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 phase 1.54, 23.16 7.30 14.85 *** .147 <.001 2 hour 1.99, 29.91 2.46 21.63 *** .099 <.001 3 phase:hour 4.10, 61.56 2.00 1.35 .011 .260 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > ### changing defaults of ANOVA table: > > # no df-correction & partial eta-squared: > aov_car(value ~ treatment * gender + Error(id/(phase*hour)), + data = obk.long, anova_table = list(correction = "none", es = "pes")) Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df MSE F pes p.value 1 treatment 2, 10 22.81 3.94 + .441 .055 2 gender 1, 10 22.81 3.66 + .268 .085 3 treatment:gender 2, 10 22.81 2.86 .364 .104 4 phase 2, 20 4.01 16.13 *** .617 <.001 5 treatment:phase 4, 20 4.01 4.85 ** .492 .007 6 gender:phase 2, 20 4.01 0.28 .028 .757 7 treatment:gender:phase 4, 20 4.01 0.64 .113 .642 8 hour 4, 40 1.56 16.69 *** .625 <.001 9 treatment:hour 8, 40 1.56 0.09 .018 .999 10 gender:hour 4, 40 1.56 0.45 .043 .772 11 treatment:gender:hour 8, 40 1.56 0.62 .110 .755 12 phase:hour 8, 80 1.20 1.18 .106 .322 13 treatment:phase:hour 16, 80 1.20 0.35 .065 .990 14 gender:phase:hour 8, 80 1.20 0.93 .085 .496 15 treatment:gender:phase:hour 16, 80 1.20 0.74 .128 .750 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > > # no df-correction and no MSE > aov_car(value ~ treatment * gender + Error(id/(phase*hour)), + data = obk.long,observed = "gender", + anova_table = list(correction = "none", MSE = FALSE)) Contrasts set to contr.sum for the following variables: treatment, gender Anova Table (Type 3 tests) Response: value Effect df F ges p.value 1 treatment 2, 10 3.94 + .198 .055 2 gender 1, 10 3.66 + .115 .085 3 treatment:gender 2, 10 2.86 .179 .104 4 phase 2, 20 16.13 *** .151 <.001 5 treatment:phase 4, 20 4.85 ** .097 .007 6 gender:phase 2, 20 0.28 .003 .757 7 treatment:gender:phase 4, 20 0.64 .014 .642 8 hour 4, 40 16.69 *** .125 <.001 9 treatment:hour 8, 40 0.09 .002 .999 10 gender:hour 4, 40 0.45 .004 .772 11 treatment:gender:hour 8, 40 0.62 .011 .755 12 phase:hour 8, 80 1.18 .015 .322 13 treatment:phase:hour 16, 80 0.35 .009 .990 14 gender:phase:hour 8, 80 0.93 .012 .496 15 treatment:gender:phase:hour 16, 80 0.74 .019 .750 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > > # add p-value adjustment for all effects (see Cramer et al., 2015, PB&R) > aov_ez("id", "value", obk.long, between = "treatment", + within = c("phase", "hour"), + anova_table = list(p_adjust_method = "holm")) Contrasts set to contr.sum for the following variables: treatment Anova Table (Type 3 tests, holm-adjusted) Response: value Effect df MSE F ges p.value 1 treatment 2, 13 32.04 2.91 .211 .360 2 phase 1.74, 22.64 4.07 19.29 *** .164 <.001 3 treatment:phase 3.48, 22.64 4.07 5.43 * .099 .022 4 hour 1.95, 25.41 2.87 18.44 *** .129 <.001 5 treatment:hour 3.91, 25.41 2.87 0.08 .001 >.999 6 phase:hour 4.02, 52.29 2.24 1.35 .017 .796 7 treatment:phase:hour 8.05, 52.29 2.24 0.33 .008 >.999 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > > ########################### > ## 2: Follow-up Analysis ## > ########################### > > # use data as above > data(obk.long, package = "afex") > > # 1. obtain afex_aov object: > a1 <- aov_ez("id", "value", obk.long, between = c("treatment", "gender"), + within = c("phase", "hour"), observed = "gender") Contrasts set to contr.sum for the following variables: treatment, gender > > > if (requireNamespace("ggplot2") & requireNamespace("emmeans")) { + # 1b. plot data using afex_plot function, for more see: + ## vignette("afex_plot_introduction", package = "afex") + + ## default plot uses multivariate model-based CIs + afex_plot(a1, "hour", "gender", c("treatment", "phase")) + + + a1b <- aov_ez("id", "value", obk.long, between = c("treatment", "gender"), + within = c("phase", "hour"), observed = "gender", + include_aov = TRUE) + ## you can use a univariate model and CIs if you refit the model with the aov + ## slot + afex_plot(a1b, "hour", "gender", c("treatment", "phase"), + emmeans_arg = list(model = "univariate")) + + ## in a mixed between-within designs, no error-bars might be preferrable: + afex_plot(a1, "hour", "gender", c("treatment", "phase"), error = "none") + } Warning: Panel(s) show a mixed within-between-design. Error bars do not allow comparisons across all means. Suppress error bars with: error = "none" Contrasts set to contr.sum for the following variables: treatment, gender Warning: Panel(s) show a mixed within-between-design. Error bars do not allow comparisons across all means. Suppress error bars with: error = "none" > > if (requireNamespace("emmeans")) { + library("emmeans") # package emmeans needs to be attached for follow-up tests. + + # 2. obtain reference grid object (default uses multivariate model): + r1 <- emmeans(a1, ~treatment +phase) + r1 + + # 3. create list of contrasts on the reference grid: + c1 <- list( + A_B_pre = c(rep(0, 6), 0, -1, 1), # A versus B for pretest + A_B_comb = c(-0.5, 0.5, 0, -0.5, 0.5, 0, 0, 0, 0), # A vs. B for post and follow-up combined + effect_post = c(0, 0, 0, -1, 0.5, 0.5, 0, 0, 0), # control versus A&B post + effect_fup = c(-1, 0.5, 0.5, 0, 0, 0, 0, 0, 0), # control versus A&B follow-up + effect_comb = c(-0.5, 0.25, 0.25, -0.5, 0.25, 0.25, 0, 0, 0) # control versus A&B combined + ) + + # 4. test contrasts on reference grid: + contrast(r1, c1) + + # same as before, but using Bonferroni-Holm correction for multiple testing: + contrast(r1, c1, adjust = "holm") + + # 2. (alternative): all pairwise comparisons of treatment: + emmeans(a1, "treatment", contr = "pairwise") + } NOTE: Results may be misleading due to involvement in interactions $emmeans treatment emmean SE df lower.CL upper.CL control 4.22 0.563 10 2.97 5.48 A 6.25 0.617 10 4.88 7.62 B 6.03 0.471 10 4.98 7.08 Results are averaged over the levels of: gender, hour, phase Confidence level used: 0.95 $contrasts contrast estimate SE df t.ratio p.value control - A -2.028 0.835 10 -2.429 0.0830 control - B -1.806 0.734 10 -2.461 0.0789 A - B 0.222 0.776 10 0.286 0.9560 Results are averaged over the levels of: gender, hour, phase P value adjustment: tukey method for comparing a family of 3 estimates > > ####################### > ## 3: Other examples ## > ####################### > data(obk.long, package = "afex") > > # replicating ?Anova using aov_car: > obk_anova <- aov_car(value ~ treatment * gender + Error(id/(phase*hour)), + data = obk.long, type = 2) Contrasts set to contr.sum for the following variables: treatment, gender > # in contrast to aov you do not need the within-subject factors outside Error() > > str(obk_anova, 1, give.attr = FALSE) List of 5 $ anova_table:Classes ‘anova’ and 'data.frame': 15 obs. of 6 variables: $ aov : NULL $ Anova :List of 14 $ lm :List of 13 $ data :List of 3 > # List of 5 > # $ anova_table:Classes ‘anova’ and 'data.frame': 15 obs. of 6 variables: > # $ aov : NULL > # $ Anova :List of 14 > # $ lm :List of 13 > # $ data :List of 3 > > obk_anova$Anova Type II Repeated Measures MANOVA Tests: Pillai test statistic Df test stat approx F num Df den Df Pr(>F) (Intercept) 1 0.96954 318.34 1 10 6.532e-09 *** treatment 2 0.48092 4.63 2 10 0.0376868 * gender 1 0.20356 2.56 1 10 0.1409735 treatment:gender 2 0.36350 2.86 2 10 0.1044692 phase 1 0.85052 25.61 2 9 0.0001930 *** treatment:phase 2 0.68518 2.61 4 20 0.0667354 . gender:phase 1 0.04314 0.20 2 9 0.8199968 treatment:gender:phase 2 0.31060 0.92 4 20 0.4721498 hour 1 0.93468 25.04 4 7 0.0003043 *** treatment:hour 2 0.30144 0.35 8 16 0.9295212 gender:hour 1 0.29274 0.72 4 7 0.6023742 treatment:gender:hour 2 0.57022 0.80 8 16 0.6131884 phase:hour 1 0.54958 0.46 8 3 0.8324517 treatment:phase:hour 2 0.66367 0.25 16 8 0.9914415 gender:phase:hour 1 0.69505 0.85 8 3 0.6202076 treatment:gender:phase:hour 2 0.79277 0.33 16 8 0.9723693 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > # Type II Repeated Measures MANOVA Tests: Pillai test statistic > # Df test stat approx F num Df den Df Pr(>F) > # (Intercept) 1 0.96954 318.34 1 10 6.532e-09 *** > # treatment 2 0.48092 4.63 2 10 0.0376868 * > # gender 1 0.20356 2.56 1 10 0.1409735 > # treatment:gender 2 0.36350 2.86 2 10 0.1044692 > # phase 1 0.85052 25.61 2 9 0.0001930 *** > # treatment:phase 2 0.68518 2.61 4 20 0.0667354 . > # gender:phase 1 0.04314 0.20 2 9 0.8199968 > # treatment:gender:phase 2 0.31060 0.92 4 20 0.4721498 > # hour 1 0.93468 25.04 4 7 0.0003043 *** > # treatment:hour 2 0.30144 0.35 8 16 0.9295212 > # gender:hour 1 0.29274 0.72 4 7 0.6023742 > # treatment:gender:hour 2 0.57022 0.80 8 16 0.6131884 > # phase:hour 1 0.54958 0.46 8 3 0.8324517 > # treatment:phase:hour 2 0.66367 0.25 16 8 0.9914415 > # gender:phase:hour 1 0.69505 0.85 8 3 0.6202076 > # treatment:gender:phase:hour 2 0.79277 0.33 16 8 0.9723693 > # --- > # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > > > cleanEx() detaching ‘package:emmeans’ > nameEx("compare.2.vectors") > ### * compare.2.vectors > > flush(stderr()); flush(stdout()) > > ### Name: compare.2.vectors > ### Title: Compare two vectors using various tests. > ### Aliases: compare.2.vectors > > ### ** Examples > > > with(sleep, compare.2.vectors(extra[group == 1], extra[group == 2])) Warning in wilcox.test.default(x, y, paired = paired, exact = wilcox.exact, : cannot compute exact p-value with ties $parametric test test.statistic test.value test.df p 1 t t -1.860813 18.00000 0.07918671 2 Welch t -1.860813 17.77647 0.07939414 $nonparametric test test.statistic test.value test.df p 1 stats::Wilcoxon W 25.500000 NA 0.06932758 2 permutation Z -1.750807 NA 0.08045000 3 coin::Wilcoxon Z -1.854118 NA 0.06665000 4 median Z -1.743560 NA 0.17750000 > > # gives: > ## $parametric > ## test test.statistic test.value test.df p > ## 1 t t -1.861 18.00 0.07919 > ## 2 Welch t -1.861 17.78 0.07939 > ## > ## $nonparametric > ## test test.statistic test.value test.df p > ## 1 stats::Wilcoxon W 25.500 NA 0.06933 > ## 2 permutation Z -1.751 NA 0.08154 > ## 3 coin::Wilcoxon Z -1.854 NA 0.06487 > ## 4 median Z -1.744 NA 0.17867 > > # compare with: > with(sleep, compare.2.vectors(extra[group == 1], extra[group == 2], + alternative = "less")) Warning in wilcox.test.default(x, y, paired = paired, exact = wilcox.exact, : cannot compute exact p-value with ties $parametric test test.statistic test.value test.df p 1 t t -1.860813 18.00000 0.03959336 2 Welch t -1.860813 17.77647 0.03969707 $nonparametric test test.statistic test.value test.df p 1 stats::Wilcoxon W 25.500000 NA 0.03466379 2 permutation Z -1.750807 NA 0.04110000 3 coin::Wilcoxon Z -1.854118 NA 0.03207000 4 median Z -1.743560 NA 0.08964000 > > with(sleep, compare.2.vectors(extra[group == 1], extra[group == 2], + alternative = "greater")) Warning in wilcox.test.default(x, y, paired = paired, exact = wilcox.exact, : cannot compute exact p-value with ties $parametric test test.statistic test.value test.df p 1 t t -1.860813 18.00000 0.9604066 2 Welch t -1.860813 17.77647 0.9603029 $nonparametric test test.statistic test.value test.df p 1 stats::Wilcoxon W 25.500000 NA 0.9707517 2 permutation Z -1.750807 NA 0.9617000 3 coin::Wilcoxon Z -1.854118 NA 0.9708100 4 median Z -1.743560 NA 0.9881000 > > # doesn't make much sense as the data is not paired, but whatever: > with(sleep, compare.2.vectors(extra[group == 1], extra[group == 2], + paired = TRUE)) Warning in wilcox.test.default(x, y, paired = paired, exact = wilcox.exact, : cannot compute exact p-value with ties Warning in wilcox.test.default(x, y, paired = paired, exact = wilcox.exact, : cannot compute exact p-value with zeroes $parametric test test.statistic test.value test.df p 1 t t -4.062128 9 0.00283289 $nonparametric test test.statistic test.value test.df p 1 stats::Wilcoxon V 0.000000 NA 0.009090698 2 permutation Z -2.543759 NA 0.004100000 3 coin::Wilcoxon Z -2.516961 NA 0.004090000 4 median Z -2.000000 NA 0.126140000 > > # from ?t.test: > compare.2.vectors(1:10,y=c(7:20, 200)) Warning in wilcox.test.default(x, y, paired = paired, exact = wilcox.exact, : cannot compute exact p-value with ties $parametric test test.statistic test.value test.df p 1 t t -1.325921 23.0000 0.1978842 2 Welch t -1.632903 14.1646 0.1245135 $nonparametric test test.statistic test.value test.df p 1 stats::Wilcoxon W 8.000000 NA 0.0002228503 2 permutation Z -1.305464 NA 0.0998500000 3 coin::Wilcoxon Z -3.719353 NA 0.0000300000 4 median Z -3.545621 NA 0.0004700000 > > > > > cleanEx() > nameEx("ems") > ### * ems > > flush(stderr()); flush(stdout()) > > ### Name: ems > ### Title: Expected values of mean squares for factorial designs Implements > ### the Cornfield-Tukey algorithm for deriving the expected values of the > ### mean squares for factorial designs. > ### Aliases: ems > > ### ** Examples > > # 2x2 mixed anova > # A varies between-subjects, B varies within-subjects > ems(r ~ A*B*S, nested="A/S", random="S") VarianceComponent Effect e B:S A:B S B A A 1 rb rbs B 1 r ras S 1 rb A:B 1 r rs B:S 1 r e 1 attr(,"terms") [1] r a b s > > # Clark (1973) example > # random Subjects, random Words, fixed Treatments > ems(r ~ S*W*T, nested="T/W", random="SW") VarianceComponent Effect e S:T S:W T W S S 1 r rwt W 1 r rs T 1 rw r rsw rs S:W 1 r S:T 1 rw r e 1 attr(,"terms") [1] r s w t > > # EMSs for Clark design if Words are fixed > ems(r ~ S*W*T, nested="T/W", random="S") VarianceComponent Effect e S:T S:W T W S S 1 rwt W 1 r rs T 1 rw rsw S:W 1 r S:T 1 rw e 1 attr(,"terms") [1] r s w t > > > > > cleanEx() > nameEx("fhch2010") > ### * fhch2010 > > flush(stderr()); flush(stdout()) > > ### Name: fhch2010 > ### Title: Data from Freeman, Heathcote, Chalmers, & Hockley (2010) > ### Aliases: fhch2010 > ### Keywords: dataset > > ### ** Examples > > > data("fhch2010") > str(fhch2010) 'data.frame': 13222 obs. of 10 variables: $ id : Factor w/ 45 levels "N1","N12","N13",..: 1 1 1 1 1 1 1 1 1 1 ... $ task : Factor w/ 2 levels "naming","lexdec": 1 1 1 1 1 1 1 1 1 1 ... $ stimulus : Factor w/ 2 levels "word","nonword": 1 1 1 2 2 1 2 2 1 2 ... $ density : Factor w/ 2 levels "low","high": 2 1 1 2 1 2 1 1 1 1 ... $ frequency: Factor w/ 2 levels "low","high": 1 2 2 2 2 2 1 2 1 2 ... $ length : Factor w/ 3 levels "4","5","6": 3 3 2 2 1 1 3 2 1 3 ... $ item : Factor w/ 600 levels "abide","acts",..: 363 121 202 525 580 135 42 368 227 141 ... $ rt : num 1.091 0.876 0.71 1.21 0.843 ... $ log_rt : num 0.0871 -0.1324 -0.3425 0.1906 -0.1708 ... $ correct : logi TRUE TRUE TRUE TRUE TRUE TRUE ... > > a1 <- aov_ez("id", "log_rt", fhch2010, between = "task", + within = c("density", "frequency", "length", "stimulus")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: task > nice(a1) Anova Table (Type 3 tests) Response: log_rt Effect df MSE F ges 1 task 1, 43 0.92 13.68 *** .202 2 density 1, 43 0.01 0.35 <.001 3 task:density 1, 43 0.01 47.23 *** .005 4 frequency 1, 43 0.02 0.84 <.001 5 task:frequency 1, 43 0.02 105.47 *** .034 6 length 1.92, 82.45 0.01 20.11 *** .007 7 task:length 1.92, 82.45 0.01 1.17 <.001 8 stimulus 1, 43 0.06 154.78 *** .155 9 task:stimulus 1, 43 0.06 78.94 *** .086 10 density:frequency 1, 43 0.00 0.47 <.001 11 task:density:frequency 1, 43 0.00 19.27 *** .002 12 density:length 1.99, 85.54 0.01 1.53 <.001 13 task:density:length 1.99, 85.54 0.01 2.41 + <.001 14 frequency:length 1.85, 79.64 0.01 2.55 + <.001 15 task:frequency:length 1.85, 79.64 0.01 3.85 * .001 16 density:stimulus 1, 43 0.01 0.40 <.001 17 task:density:stimulus 1, 43 0.01 25.71 *** .003 18 frequency:stimulus 1, 43 0.01 114.85 *** .017 19 task:frequency:stimulus 1, 43 0.01 192.62 *** .028 20 length:stimulus 1.71, 73.58 0.01 1.42 <.001 21 task:length:stimulus 1.71, 73.58 0.01 1.07 <.001 22 density:frequency:length 1.99, 85.68 0.01 1.30 <.001 23 task:density:frequency:length 1.99, 85.68 0.01 1.94 <.001 24 density:frequency:stimulus 1, 43 0.01 3.54 + <.001 25 task:density:frequency:stimulus 1, 43 0.01 10.33 ** .002 26 density:length:stimulus 1.84, 79.10 0.01 0.54 <.001 27 task:density:length:stimulus 1.84, 79.10 0.01 0.36 <.001 28 frequency:length:stimulus 1.90, 81.51 0.01 0.40 <.001 29 task:frequency:length:stimulus 1.90, 81.51 0.01 11.06 *** .003 30 density:frequency:length:stimulus 1.86, 79.95 0.01 0.87 <.001 31 task:density:frequency:length:stimulus 1.86, 79.95 0.01 0.00 <.001 p.value 1 <.001 2 .559 3 <.001 4 .364 5 <.001 6 <.001 7 .314 8 <.001 9 <.001 10 .495 11 <.001 12 .223 13 .096 14 .088 15 .028 16 .529 17 <.001 18 <.001 19 <.001 20 .248 21 .340 22 .278 23 .151 24 .067 25 .002 26 .571 27 .680 28 .658 29 <.001 30 .416 31 .994 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > if (requireNamespace("emmeans") && requireNamespace("ggplot2")) { + afex_plot(a1, "length", "frequency", c("task", "stimulus"), error = "within") + + afex_plot(a1, "density", "frequency", c("task", "stimulus"), error = "within") + } > > > ## Not run: > ##D a2 <- aov_ez("id", "rt", fhch2010, between = "task", > ##D within = c("density", "frequency", "length", "stimulus")) > ##D nice(a2) > ##D > ##D if (requireNamespace("emmeans") && requireNamespace("ggplot2")) { > ##D afex_plot(a2, "length", "frequency", c("task", "stimulus"), error = "within") > ##D > ##D afex_plot(a2, "density", "frequency", c("task", "stimulus"), error = "within") > ##D } > ## End(Not run) > > > > cleanEx() > nameEx("ks2013.3") > ### * ks2013.3 > > flush(stderr()); flush(stdout()) > > ### Name: ks2013.3 > ### Title: Data from Klauer & Singmann (2013, Experiment 3) > ### Aliases: ks2013.3 > ### Keywords: dataset > > ### ** Examples > > data("ks2013.3") > > # replicate results reported in Klauer & Singmann (2013, p. 1270) > > aov_ez("id", "response", ks2013.3, between = "condition", + within = c("believability", "validity")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: condition Anova Table (Type 3 tests) Response: response Effect df MSE F ges p.value 1 condition 1, 58 0.94 0.01 <.001 .903 2 believability 1.84, 106.78 0.59 8.36 *** .051 <.001 3 condition:believability 1.84, 106.78 0.59 0.29 .002 .728 4 validity 1, 58 0.38 0.17 <.001 .684 5 condition:validity 1, 58 0.38 2.07 .005 .155 6 believability:validity 1.85, 107.52 0.28 8.29 *** .025 <.001 7 condition:believability:validity 1.85, 107.52 0.28 3.58 * .011 .035 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_ez("id", "response", subset(ks2013.3, condition == "fixed"), + within = c("believability", "validity")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Anova Table (Type 3 tests) Response: response Effect df MSE F ges p.value 1 believability 1.58, 45.92 0.79 4.61 * .057 .022 2 validity 1, 29 0.31 2.10 .007 .158 3 believability:validity 1.75, 50.88 0.34 8.66 *** .051 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_ez("id", "response", subset(ks2013.3, condition == "random"), + within = c("believability", "validity")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Anova Table (Type 3 tests) Response: response Effect df MSE F ges p.value 1 believability 1.87, 54.28 0.48 3.94 * .046 .028 2 validity 1, 29 0.45 0.45 .003 .508 3 believability:validity 1.94, 56.31 0.23 2.25 .013 .116 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > > > cleanEx() > nameEx("md_12.1") > ### * md_12.1 > > flush(stderr()); flush(stdout()) > > ### Name: md_12.1 > ### Title: Data 12.1 from Maxwell & Delaney > ### Aliases: md_12.1 > ### Keywords: dataset > > ### ** Examples > > data(md_12.1) > > # Table 12.5 (p. 578): > aov_ez("id", "rt", md_12.1, within = c("angle", "noise"), + args.return=list(correction = "none", es = "none")) Anova Table (Type 3 tests) Response: rt Effect df MSE F ges p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .390 <.001 2 noise 1, 9 8460.00 33.77 *** .387 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .188 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > > > > > cleanEx() > nameEx("md_15.1") > ### * md_15.1 > > flush(stderr()); flush(stdout()) > > ### Name: md_15.1 > ### Title: Data 15.1 / 11.5 from Maxwell & Delaney > ### Aliases: md_15.1 > ### Keywords: dataset > > ### ** Examples > > ### replicate results from Table 15.2 to 15.6 (Maxwell & Delaney, 2004, pp. 774) > data(md_15.1) > > ### ANOVA results (Table 15.2) > aov_4(iq ~ timecat + (timecat|id),data=md_15.1, anova_table=list(correction = "none")) Anova Table (Type 3 tests) Response: iq Effect df MSE F ges p.value 1 timecat 3, 33 60.79 3.03 * .060 .043 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > > ### Table 15.3 (random intercept only) > # we need to set the base level on the last level: > contrasts(md_15.1$timecat) <- contr.treatment(4, base = 4) > # "Type 3 Tests of Fixed Effects" > (t15.3 <- mixed(iq ~ timecat + (1|id),data=md_15.1, check.contrasts=FALSE)) Warning: 'check.contrasts' is deprecated; use 'check_contrasts' instead Mixed Model Anova Table (Type 3 tests, S-method) Model: iq ~ timecat + (1 | id) Data: md_15.1 Effect df F p.value 1 timecat 3, 33.00 3.03 * .043 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > # "Solution for Fixed Effects" and "Covariance Parameter Estimates" > summary(t15.3$full.model) Length Class Mode 0 NULL NULL > > ### make Figure 15.2 > plot(NULL, NULL, ylim = c(80, 140), xlim = c(30, 48), ylab = "iq", xlab = "time") > plyr::d_ply(md_15.1, plyr::.(id), function(x) lines(as.numeric(as.character(x$timecat)), x$iq)) > > ### Table 15.4, page 789 > # random intercept plus slope > (t15.4 <- mixed(iq ~ timecat + (1+time|id),data=md_15.1, check.contrasts=FALSE)) Warning: 'check.contrasts' is deprecated; use 'check_contrasts' instead Mixed Model Anova Table (Type 3 tests, S-method) Model: iq ~ timecat + (1 + time | id) Data: md_15.1 Effect df F p.value 1 timecat 3, 16.73 1.78 .189 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > summary(t15.4$full.model) Length Class Mode 0 NULL NULL > > ### Table 15.5, page 795 > # set up polynomial contrasts for timecat > contrasts(md_15.1$timecat) <- contr.poly > # fit all parameters separately > (t15.5 <- mixed(iq ~ timecat + (1+time|id), data=md_15.1, check.contrasts=FALSE, + per.parameter="timecat")) Warning: 'per.parameter' is deprecated; use 'per_parameter' instead Warning: 'check.contrasts' is deprecated; use 'check_contrasts' instead Mixed Model Anova Table (Type 3 tests, S-method) Model: iq ~ timecat + (1 + time | id) Data: md_15.1 Effect df F p.value 1 timecat 3, 16.21 1.78 .191 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > # quadratic trend is considerably off, conclusions stay the same. > > > ### Table 15.6, page 797 > # growth curve model > (t15.6 <- mixed(iq ~ time + (1+time|id),data=md_15.1)) Contrasts set to contr.sum for the following variables: id Numerical variables NOT centered on 0: time If in interactions, interpretation of lower order (e.g., main) effects difficult. Mixed Model Anova Table (Type 3 tests, S-method) Model: iq ~ time + (1 + time | id) Data: md_15.1 Effect df F p.value 1 time 1, 11.00 5.02 * .047 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > summary(t15.6$full.model) Length Class Mode 0 NULL NULL > > > > > cleanEx() > nameEx("md_16.1") > ### * md_16.1 > > flush(stderr()); flush(stdout()) > > ### Name: md_16.1 > ### Title: Data 16.1 / 10.9 from Maxwell & Delaney > ### Aliases: md_16.1 > ### Keywords: dataset > > ### ** Examples > > ### replicate results from Table 16.3 (Maxwell & Delaney, 2004, p. 837) > data(md_16.1) > > # original results need treatment contrasts: > (mixed1_orig <- mixed(severity ~ sex + (1|id), md_16.1, check.contrasts=FALSE)) Warning: 'check.contrasts' is deprecated; use 'check_contrasts' instead Mixed Model Anova Table (Type 3 tests, S-method) Model: severity ~ sex + (1 | id) Data: md_16.1 Effect df F p.value 1 sex 1, 4.00 9.97 * .034 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > summary(mixed1_orig$full.model) Length Class Mode 0 NULL NULL > > # p-values stay the same with afex default contrasts (contr.sum), > # but estimates and t-values for the fixed effects parameters change. > (mixed1 <- mixed(severity ~ sex + (1|id), md_16.1)) Contrasts set to contr.sum for the following variables: sex, id Mixed Model Anova Table (Type 3 tests, S-method) Model: severity ~ sex + (1 | id) Data: md_16.1 Effect df F p.value 1 sex 1, 4.00 9.97 * .034 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > summary(mixed1$full.model) Length Class Mode 0 NULL NULL > > > > > cleanEx() > nameEx("md_16.4") > ### * md_16.4 > > flush(stderr()); flush(stdout()) > > ### Name: md_16.4 > ### Title: Data 16.4 from Maxwell & Delaney > ### Aliases: md_16.4 > ### Keywords: dataset > > ### ** Examples > > # data for next examples (Maxwell & Delaney, Table 16.4) > data(md_16.4) > str(md_16.4) 'data.frame': 29 obs. of 6 variables: $ obs : Factor w/ 29 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ... $ room : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 2 2 ... $ cond : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ... $ cog : int 46 52 60 44 46 48 50 54 50 48 ... $ skill : int 4 4 4 4 6 6 6 6 6 6 ... $ induct: int 21 26 33 22 18 25 26 24 21 25 ... > > ### replicate results from Table 16.6 (Maxwell & Delaney, 2004, p. 845) > # p-values (almost) hold: > (mixed2 <- mixed(induct ~ cond + (1|room:cond), md_16.4)) Contrasts set to contr.sum for the following variables: cond, room Mixed Model Anova Table (Type 3 tests, S-method) Model: induct ~ cond + (1 | room:cond) Data: md_16.4 Effect df F p.value 1 cond 1, 3.90 3.21 .150 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > # (1|room:cond) is needed because room is nested within cond. > > > > > > cleanEx() > nameEx("mixed") > ### * mixed > > flush(stderr()); flush(stdout()) > > ### Encoding: UTF-8 > > ### Name: mixed > ### Title: p-values for fixed effects of mixed-model via lme4::lmer() > ### Aliases: mixed lmer_alt > > ### ** Examples > > > ################################## > ## Simple Examples (from MEMSS) ## > ################################## > > if (requireNamespace("MEMSS")) { + data("Machines", package = "MEMSS") + + # simple model with random-slopes for repeated-measures factor + m1 <- mixed(score ~ Machine + (Machine|Worker), data=Machines) + m1 + + # suppress correlations among random effect parameters with || and expand_re = TRUE + m2 <- mixed(score ~ Machine + (Machine||Worker), data=Machines, expand_re = TRUE) + m2 + + ## compare: + summary(m1)$varcor + summary(m2)$varcor + # for wrong solution see: + # summary(lmer(score ~ Machine + (Machine||Worker), data=Machines))$varcor + + if (requireNamespace("emmeans")) { + # follow-up tests + library("emmeans") # package emmeans needs to be attached for follow-up tests. + (emm1 <- emmeans(m1, "Machine")) + pairs(emm1, adjust = "holm") # all pairwise comparisons + con1 <- list( + c1 = c(1, -0.5, -0.5), # 1 versus other 2 + c2 = c(0.5, -1, 0.5) # 1 and 3 versus 2 + ) + contrast(emm1, con1, adjust = "holm") + + if (requireNamespace("ggplot2")) { + # plotting + afex_plot(m1, "Machine") ## default uses model-based CIs + ## within-subjects CIs somewhat more in line with pairwirse comparisons: + afex_plot(m1, "Machine", error = "within") + + ## less differences between CIs for model without correlations: + afex_plot(m2, "Machine") + afex_plot(m2, "Machine", error = "within") + }}} Contrasts set to contr.sum for the following variables: Machine, Worker Contrasts set to contr.sum for the following variables: Machine, Worker Aggregating data over: Worker Aggregating data over: Worker Aggregating data over: Worker Aggregating data over: Worker > > ## Not run: > ##D ####################### > ##D ### Further Options ### > ##D ####################### > ##D > ##D ## Multicore: > ##D > ##D require(parallel) > ##D (nc <- detectCores()) # number of cores > ##D cl <- makeCluster(rep("localhost", nc)) # make cluster > ##D # to keep track of what the function is doindg redirect output to outfile: > ##D # cl <- makeCluster(rep("localhost", nc), outfile = "cl.log.txt") > ##D > ##D data("Machines", package = "MEMSS") > ##D ## There are two ways to use multicore: > ##D > ##D # 1. Obtain fits with multicore (e.g. for likelihood ratio tests, LRT): > ##D mixed(score ~ Machine + (Machine|Worker), data=Machines, cl = cl, > ##D method = "LRT") > ##D > ##D # 2. Obtain PB samples via multicore: > ##D mixed(score ~ Machine + (Machine|Worker), data=Machines, > ##D method = "PB", args_test = list(nsim = 50, cl = cl)) # better use 500 or 1000 > ##D > ##D ## Both ways can be combined: > ##D # 2. Obtain PB samples via multicore: > ##D mixed(score ~ Machine + (Machine|Worker), data=Machines, cl = cl, > ##D method = "PB", args_test = list(nsim = 50, cl = cl)) > ##D > ##D #### use all_fit = TRUE and expand_re = TRUE: > ##D data("sk2011.2") # data described in more detail below > ##D sk2_aff <- droplevels(sk2011.2[sk2011.2$what == "affirmation",]) > ##D > ##D require(optimx) # uses two more algorithms > ##D sk2_aff_b <- mixed(response ~ instruction*type+(inference*type||id), sk2_aff, > ##D expand_re = TRUE, all_fit = TRUE, method = "LRT") > ##D attr(sk2_aff_b, "all_fit_selected") > ##D attr(sk2_aff_b, "all_fit_logLik") > ##D > ##D # considerably faster with multicore: > ##D clusterEvalQ(cl, library(optimx)) # need to load optimx in cluster > ##D sk2_aff_b2 <- mixed(response ~ instruction*type+(inference*type||id), sk2_aff, > ##D expand_re = TRUE, all_fit = TRUE, cl=cl, method = "LRT") > ##D attr(sk2_aff_b2, "all_fit_selected") > ##D attr(sk2_aff_b2, "all_fit_logLik") > ##D > ##D > ##D stopCluster(cl) > ##D > ## End(Not run) > > ################################################### > ## Replicating Maxwell & Delaney (2004) Examples ## > ################################################### > ## Not run: > ##D > ##D ### replicate results from Table 15.4 (Maxwell & Delaney, 2004, p. 789) > ##D data(md_15.1) > ##D # random intercept plus random slope > ##D (t15.4a <- mixed(iq ~ timecat + (1+time|id),data=md_15.1)) > ##D > ##D # to also replicate exact parameters use treatment.contrasts and the last level as base level: > ##D contrasts(md_15.1$timecat) <- contr.treatment(4, base = 4) > ##D (t15.4b <- mixed(iq ~ timecat + (1+time|id),data=md_15.1, check_contrasts=FALSE)) > ##D summary(t15.4a) # gives "wrong" parameters extimates > ##D summary(t15.4b) # identical parameters estimates > ##D > ##D # for more examples from chapter 15 see ?md_15.1 > ##D > ##D ### replicate results from Table 16.3 (Maxwell & Delaney, 2004, p. 837) > ##D data(md_16.1) > ##D > ##D # original results need treatment contrasts: > ##D (mixed1_orig <- mixed(severity ~ sex + (1|id), md_16.1, check_contrasts=FALSE)) > ##D summary(mixed1_orig$full_model) > ##D > ##D # p-value stays the same with afex default contrasts (contr.sum), > ##D # but estimates and t-values for the fixed effects parameters change. > ##D (mixed1 <- mixed(severity ~ sex + (1|id), md_16.1)) > ##D summary(mixed1$full_model) > ##D > ##D > ##D # data for next examples (Maxwell & Delaney, Table 16.4) > ##D data(md_16.4) > ##D str(md_16.4) > ##D > ##D ### replicate results from Table 16.6 (Maxwell & Delaney, 2004, p. 845) > ##D # Note that (1|room:cond) is needed because room is nested within cond. > ##D # p-value (almost) holds. > ##D (mixed2 <- mixed(induct ~ cond + (1|room:cond), md_16.4)) > ##D # (differences are dut to the use of Kenward-Roger approximation here, > ##D # whereas M&W's p-values are based on uncorrected df.) > ##D > ##D # again, to obtain identical parameter and t-values, use treatment contrasts: > ##D summary(mixed2) # not identical > ##D > ##D # prepare new data.frame with contrasts: > ##D md_16.4b <- within(md_16.4, cond <- C(cond, contr.treatment, base = 2)) > ##D str(md_16.4b) > ##D > ##D # p-value stays identical: > ##D (mixed2_orig <- mixed(induct ~ cond + (1|room:cond), md_16.4b, > ##D check_contrasts=FALSE)) > ##D summary(mixed2_orig$full_model) # replicates parameters > ##D > ##D > ##D ### replicate results from Table 16.7 (Maxwell & Delaney, 2004, p. 851) > ##D # F-values (almost) hold, p-values (especially for skill) are off > ##D (mixed3 <- mixed(induct ~ cond + skill + (1|room:cond), md_16.4)) > ##D > ##D # however, parameters are perfectly recovered when using the original contrasts: > ##D mixed3_orig <- mixed(induct ~ cond + skill + (1|room:cond), md_16.4b, > ##D check_contrasts=FALSE) > ##D summary(mixed3_orig) > ##D > ##D > ##D ### replicate results from Table 16.10 (Maxwell & Delaney, 2004, p. 862) > ##D # for this we need to center cog: > ##D md_16.4b$cog <- scale(md_16.4b$cog, scale=FALSE) > ##D > ##D # F-values and p-values are relatively off: > ##D (mixed4 <- mixed(induct ~ cond*cog + (cog|room:cond), md_16.4b)) > ##D # contrast has a relatively important influence on cog > ##D (mixed4_orig <- mixed(induct ~ cond*cog + (cog|room:cond), md_16.4b, > ##D check_contrasts=FALSE)) > ##D > ##D # parameters are again almost perfectly recovered: > ##D summary(mixed4_orig) > ## End(Not run) > > ########################### > ## Full Analysis Example ## > ########################### > > ## Not run: > ##D ### split-plot experiment (Singmann & Klauer, 2011, Exp. 2) > ##D ## between-factor: instruction > ##D ## within-factor: inference & type > ##D ## hypothesis: three-way interaction > ##D data("sk2011.2") > ##D > ##D # use only affirmation problems (S&K also splitted the data like this) > ##D sk2_aff <- droplevels(sk2011.2[sk2011.2$what == "affirmation",]) > ##D > ##D # set up model with maximal by-participant random slopes > ##D sk_m1 <- mixed(response ~ instruction*inference*type+(inference*type|id), sk2_aff) > ##D > ##D sk_m1 # prints ANOVA table with nicely rounded numbers (i.e., as characters) > ##D nice(sk_m1) # returns the same but without printing potential warnings > ##D anova(sk_m1) # returns and prints numeric ANOVA table (i.e., not-rounded) > ##D summary(sk_m1) # lmer summary of full model > ##D > ##D # same model but using Kenward-Roger approximation of df > ##D # very similar results but slower > ##D sk_m1b <- mixed(response ~ instruction*inference*type+(inference*type|id), > ##D sk2_aff, method="KR") > ##D nice(sk_m1b) > ##D # identical results as: > ##D anova(sk_m1$full_model) > ##D > ##D # suppressing correlation among random slopes: very similar results, but > ##D # significantly faster and often less convergence warnings. > ##D sk_m2 <- mixed(response ~ instruction*inference*type+(inference*type||id), sk2_aff, > ##D expand_re = TRUE) > ##D sk_m2 > ##D > ##D ## mixed objects can be passed to emmeans > ##D library("emmeans") # however, package emmeans needs to be attached first > ##D > ##D # emmeans also approximate df which takes time with default Kenward-Roger > ##D emm_options(lmer.df = "Kenward-Roger") # default setting, slow > ##D emm_options(lmer.df = "Satterthwaite") # faster setting, preferrable > ##D emm_options(lmer.df = "asymptotic") # the fastest, df = infinity > ##D > ##D > ##D # recreates basically Figure 4 (S&K, 2011, upper panel) > ##D # only the 4th and 6th x-axis position are flipped > ##D afex_plot(sk_m1, x = c("type", "inference"), trace = "instruction") > ##D > ##D # set up reference grid for custom contrasts: > ##D (rg1 <- emmeans(sk_m1, c("instruction", "type", "inference"))) > ##D > ##D # set up contrasts on reference grid: > ##D contr_sk2 <- list( > ##D ded_validity_effect = c(rep(0, 4), 1, rep(0, 5), -1, 0), > ##D ind_validity_effect = c(rep(0, 5), 1, rep(0, 5), -1), > ##D counter_MP = c(rep(0, 4), 1, -1, rep(0, 6)), > ##D counter_AC = c(rep(0, 10), 1, -1) > ##D ) > ##D > ##D # test the main double dissociation (see S&K, p. 268) > ##D contrast(rg1, contr_sk2, adjust = "holm") > ##D # all effects are significant. > ## End(Not run) > > #################### > ## Other Examples ## > #################### > > ## Not run: > ##D > ##D # use the obk.long data (not reasonable, no random slopes) > ##D data(obk.long) > ##D mixed(value ~ treatment * phase + (1|id), obk.long) > ##D > ##D # Examples for using the per.parameter argument > ##D # note, require method = "nested-KR", "LRT", or "PB" > ##D # also we use custom contrasts > ##D data(obk.long, package = "afex") > ##D obk.long$hour <- ordered(obk.long$hour) > ##D contrasts(obk.long$phase) <- "contr.sum" > ##D contrasts(obk.long$treatment) <- "contr.sum" > ##D > ##D # tests only the main effect parameters of hour individually per parameter. > ##D mixed(value ~ treatment*phase*hour +(1|id), per_parameter = "^hour$", > ##D data = obk.long, method = "nested-KR", check_contrasts = FALSE) > ##D > ##D # tests all parameters including hour individually > ##D mixed(value ~ treatment*phase*hour +(1|id), per_parameter = "hour", > ##D data = obk.long, method = "nested-KR", check_contrasts = FALSE) > ##D > ##D # tests all parameters individually > ##D mixed(value ~ treatment*phase*hour +(1|id), per_parameter = ".", > ##D data = obk.long, method = "nested-KR", check_contrasts = FALSE) > ##D > ##D # example data from package languageR: Lexical decision latencies elicited from > ##D # 21 subjects for 79 English concrete nouns, with variables linked to subject or > ##D # word. > ##D data(lexdec, package = "languageR") > ##D > ##D # using the simplest model > ##D m1 <- mixed(RT ~ Correct + Trial + PrevType * meanWeight + > ##D Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), data = lexdec) > ##D m1 > ##D # Mixed Model Anova Table (Type 3 tests, S-method) > ##D # > ##D # Model: RT ~ Correct + Trial + PrevType * meanWeight + Frequency + NativeLanguage * > ##D # Model: Length + (1 | Subject) + (1 | Word) > ##D # Data: lexdec > ##D # Effect df F p.value > ##D # 1 Correct 1, 1627.67 8.16 ** .004 > ##D # 2 Trial 1, 1591.92 7.58 ** .006 > ##D # 3 PrevType 1, 1605.05 0.17 .680 > ##D # 4 meanWeight 1, 74.37 14.85 *** <.001 > ##D # 5 Frequency 1, 75.06 56.54 *** <.001 > ##D # 6 NativeLanguage 1, 27.12 0.70 .412 > ##D # 7 Length 1, 74.80 8.70 ** .004 > ##D # 8 PrevType:meanWeight 1, 1600.79 6.19 * .013 > ##D # 9 NativeLanguage:Length 1, 1554.49 14.24 *** <.001 > ##D # --- > ##D # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > ##D > ##D # Fitting a GLMM using parametric bootstrap: > ##D require("mlmRev") # for the data, see ?Contraception > ##D > ##D gm1 <- mixed(use ~ age + I(age^2) + urban + livch + (1 | district), method = "PB", > ##D family = binomial, data = Contraception, args_test = list(nsim = 10)) > ##D ## note that nsim = 10 is way too low for all real examples! > ##D > ## End(Not run) > > ## Not run: > ##D ##################################### > ##D ## Interplay with effects packages ## > ##D ##################################### > ##D > ##D data("Machines", package = "MEMSS") > ##D # simple model with random-slopes for repeated-measures factor > ##D m1 <- mixed(score ~ Machine + (Machine|Worker), data=Machines, > ##D set_data_arg = TRUE) ## necessary for it to work! > ##D > ##D library("effects") > ##D > ##D Effect("Machine", m1$full_model) # not correct: > ##D # Machine effect > ##D # Machine > ##D # A B C > ##D # 59.65000 52.35556 60.32222 > ##D > ##D # compare: > ##D emmeans::emmeans(m1, "Machine") > ##D # Machine emmean SE df asymp.LCL asymp.UCL > ##D # A 52.35556 1.680711 Inf 49.06142 55.64969 > ##D # B 60.32222 3.528546 Inf 53.40640 67.23804 > ##D # C 66.27222 1.806273 Inf 62.73199 69.81245 > ##D > ##D ## necessary to set contr.sum globally: > ##D set_sum_contrasts() > ##D Effect("Machine", m1$full_model) > ##D # Machine effect > ##D # Machine > ##D # A B C > ##D # 52.35556 60.32222 66.27222 > ##D > ##D plot(Effect("Machine", m1$full_model)) > ## End(Not run) > > > > cleanEx() detaching ‘package:emmeans’ > nameEx("nice") > ### * nice > > flush(stderr()); flush(stdout()) > > ### Name: nice > ### Title: Make nice ANOVA table for printing. > ### Aliases: nice nice.afex_aov nice.anova nice.mixed print.nice_table > > ### ** Examples > > > ## example from Olejnik & Algina (2003) > # "Repeated Measures Design" (pp. 439): > data(md_12.1) > # create object of class afex_aov: > rmd <- aov_ez("id", "rt", md_12.1, within = c("angle", "noise")) > rmd Anova Table (Type 3 tests) Response: rt Effect df MSE F ges p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .390 <.001 2 noise 1, 9 8460.00 33.77 *** .387 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .188 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > nice(rmd) Anova Table (Type 3 tests) Response: rt Effect df MSE F ges p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .390 <.001 2 noise 1, 9 8460.00 33.77 *** .387 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .188 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > str(nice(rmd)) Classes ‘nice_table’ and 'data.frame': 3 obs. of 6 variables: $ Effect : chr "angle" "noise" "angle:noise" $ df : chr "1.92, 17.31" "1, 9" "1.81, 16.27" $ MSE : chr "3702.02" "8460.00" "1283.22" $ F : chr "40.72 ***" "33.77 ***" "45.31 ***" $ ges : chr ".390" ".387" ".188" $ p.value: chr "<.001" "<.001" "<.001" - attr(*, "heading")= chr [1:2] "Anova Table (Type 3 tests)\n" "Response: rt" - attr(*, "p_adjust_method")= chr "none" - attr(*, "correction")= chr "GG" - attr(*, "observed")= chr(0) - attr(*, "es")= chr "ges" - attr(*, "sig_symbols")= chr [1:4] " +" " *" " **" " ***" - attr(*, "round_ps")=function (x) > # use different es: > nice(rmd, es = "pes") # noise: .82 Anova Table (Type 3 tests) Response: rt Effect df MSE F pes p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .819 <.001 2 noise 1, 9 8460.00 33.77 *** .790 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .834 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > nice(rmd, es = "ges") # noise: .39 Anova Table (Type 3 tests) Response: rt Effect df MSE F ges p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .390 <.001 2 noise 1, 9 8460.00 33.77 *** .387 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .188 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > # same data other approach: > rmd2 <- aov_ez("id", "rt", md_12.1, within = c("angle", "noise"), + anova_table=list(correction = "none", es = "none")) > nice(rmd2) Anova Table (Type 3 tests) Response: rt Effect df MSE F p.value 1 angle 2, 18 3560.00 40.72 *** <.001 2 noise 1, 9 8460.00 33.77 *** <.001 3 angle:noise 2, 18 1160.00 45.31 *** <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > nice(rmd2, correction = "GG") Anova Table (Type 3 tests) Response: rt Effect df MSE F p.value 1 angle 1.92, 17.31 3702.02 40.72 *** <.001 2 noise 1, 9 8460.00 33.77 *** <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > nice(rmd2, correction = "GG", es = "ges") Anova Table (Type 3 tests) Response: rt Effect df MSE F ges p.value 1 angle 1.92, 17.31 3702.02 40.72 *** .390 <.001 2 noise 1, 9 8460.00 33.77 *** .387 <.001 3 angle:noise 1.81, 16.27 1283.22 45.31 *** .188 <.001 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > # exampel using obk.long (see ?obk.long), a long version of the OBrienKaiser dataset from car. > data(obk.long) > # create object of class afex_aov: > tmp.aov <- aov_car(value ~ treatment * gender + Error(id/phase*hour), data = obk.long) Contrasts set to contr.sum for the following variables: treatment, gender > > nice(tmp.aov, observed = "gender") Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 22.81 3.94 + .198 .055 2 gender 1, 10 22.81 3.66 + .115 .085 3 treatment:gender 2, 10 22.81 2.86 .179 .104 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 8 hour 1.84, 18.41 3.39 16.69 *** .125 <.001 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > nice(tmp.aov, observed = "gender", sig_symbols = rep("", 4)) Anova Table (Type 3 tests) Response: value Effect df MSE F ges p.value 1 treatment 2, 10 22.81 3.94 .198 .055 2 gender 1, 10 22.81 3.66 .115 .085 3 treatment:gender 2, 10 22.81 2.86 .179 .104 4 phase 1.60, 15.99 5.02 16.13 .151 <.001 5 treatment:phase 3.20, 15.99 5.02 4.85 .097 .013 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612 8 hour 1.84, 18.41 3.39 16.69 .125 <.001 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646 Sphericity correction method: GG > > ## Not run: > ##D # use package ascii or xtable for formatting of tables ready for printing. > ##D > ##D full <- nice(tmp.aov, observed = "gender") > ##D > ##D require(ascii) > ##D print(ascii(full, include.rownames = FALSE, caption = "ANOVA 1"), type = "org") > ##D > ##D require(xtable) > ##D print.xtable(xtable(full, caption = "ANOVA 2"), include.rownames = FALSE) > ## End(Not run) > > > > cleanEx() > nameEx("obk.long") > ### * obk.long > > flush(stderr()); flush(stdout()) > > ### Name: obk.long > ### Title: O'Brien Kaiser's Repeated-Measures Dataset with Covariate > ### Aliases: obk.long > ### Keywords: dataset > > ### ** Examples > > # The dataset is constructed as follows: > data("OBrienKaiser", package = "carData") > set.seed(1) > OBrienKaiser2 <- within(OBrienKaiser, { + id <- factor(1:nrow(OBrienKaiser)) + age <- scale(sample(18:35, nrow(OBrienKaiser), replace = TRUE), scale = FALSE)}) > attributes(OBrienKaiser2$age) <- NULL # needed or resahpe2::melt throws an error. > OBrienKaiser2$age <- as.numeric(OBrienKaiser2$age) > obk.long <- reshape2::melt(OBrienKaiser2, id.vars = c("id", "treatment", "gender", "age")) > obk.long[,c("phase", "hour")] <- lapply(as.data.frame(do.call(rbind, + strsplit(as.character(obk.long$variable), "\\."),)), factor) > obk.long <- obk.long[,c("id", "treatment", "gender", "age", "phase", "hour", "value")] > obk.long <- obk.long[order(obk.long$id),] > rownames(obk.long) <- NULL > str(obk.long) 'data.frame': 240 obs. of 7 variables: $ id : Factor w/ 16 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ... $ treatment: Factor w/ 3 levels "control","A",..: 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "contrasts")= num [1:3, 1:2] -2 1 1 0 -1 1 .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr [1:3] "control" "A" "B" .. .. ..$ : NULL $ gender : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ... ..- attr(*, "contrasts")= chr "contr.sum" $ age : num -4.56 -4.56 -4.56 -4.56 -4.56 ... $ phase : Factor w/ 3 levels "fup","post","pre": 3 3 3 3 3 2 2 2 2 2 ... $ hour : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ... $ value : num 1 2 4 2 1 3 2 5 3 2 ... > ## 'data.frame': 240 obs. of 7 variables: > ## $ id : Factor w/ 16 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ... > ## $ treatment: Factor w/ 3 levels "control","A",..: 1 1 1 1 1 1 1 1 1 1 ... > ## $ gender : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ... > ## $ age : num -4.75 -4.75 -4.75 -4.75 -4.75 -4.75 -4.75 -4.75 -4.75 -4.75 ... > ## $ phase : Factor w/ 3 levels "fup","post","pre": 3 3 3 3 3 2 2 2 2 2 ... > ## $ hour : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ... > ## $ value : num 1 2 4 2 1 3 2 5 3 2 ... > head(obk.long) id treatment gender age phase hour value 1 1 control M -4.5625 pre 1 1 2 1 control M -4.5625 pre 2 2 3 1 control M -4.5625 pre 3 4 4 1 control M -4.5625 pre 4 2 5 1 control M -4.5625 pre 5 1 6 1 control M -4.5625 post 1 3 > ## id treatment gender age phase hour value > ## 1 1 control M -4.75 pre 1 1 > ## 2 1 control M -4.75 pre 2 2 > ## 3 1 control M -4.75 pre 3 4 > ## 4 1 control M -4.75 pre 4 2 > ## 5 1 control M -4.75 pre 5 1 > ## 6 1 control M -4.75 post 1 3 > > > > cleanEx() > nameEx("predict.afex_aov") > ### * predict.afex_aov > > flush(stderr()); flush(stdout()) > > ### Name: predict.afex_aov > ### Title: Predict method for 'afex_aov' objects > ### Aliases: predict.afex_aov > > ### ** Examples > > > data(obk.long, package = "afex") > > # estimate mixed ANOVA on the full design: > fit <- aov_ez("id", "value", obk.long, between = c("treatment", "gender"), + within = c("phase", "hour"), observed = "gender") Contrasts set to contr.sum for the following variables: treatment, gender > > new_data <- expand.grid( + treatment = "A", + gender = "F", + phase = c("pre", "post"), + hour = c(1, 5) + ) > > predict(fit, newdata = new_data) 1 2 3 4 2.5 3.0 3.0 3.0 > predict(fit, newdata = new_data, append = TRUE) treatment gender phase hour .predict 1 A F pre 1 2.5 2 A F post 1 3.0 3 A F pre 5 3.0 4 A F post 5 3.0 > > > > cleanEx() > nameEx("residuals.afex_aov") > ### * residuals.afex_aov > > flush(stderr()); flush(stdout()) > > ### Name: residuals.afex_aov > ### Title: Extract Residuals and Fitted Values from 'afex_aov' objects > ### Aliases: residuals.afex_aov fitted.afex_aov > > ### ** Examples > > ### Setup ANOVAs > data(obk.long, package = "afex") > between <- aov_car(value ~ treatment*gender + Error(id), data = obk.long) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: treatment, gender > within <- aov_car(value ~ 1 + Error(id/(phase*hour)), data = obk.long) > mixed <- aov_car(value ~ treatment * gender + Error(id/(phase*hour)), data = obk.long) Contrasts set to contr.sum for the following variables: treatment, gender > > # All residuals call produce the message that the data was changed during calculation. > residuals(within) Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data. residuals(..., append = TRUE) will return data and residuals. 1 2 3 4 5 6 7 8 9 10 -4.0625 -3.2500 -5.8125 -2.3125 -1.4375 -2.0625 -3.5625 -1.8750 -3.3750 -2.8750 11 12 13 14 15 16 17 18 19 20 -2.8125 -2.1875 -1.3125 -2.5000 -3.0625 1.9375 1.7500 1.1875 0.6875 2.5625 21 22 23 24 25 26 27 28 29 30 0.9375 1.4375 -0.8750 1.6250 3.1250 0.1875 -0.1875 -0.3125 -1.5000 -0.0625 31 32 33 34 35 36 37 38 39 40 -1.0625 -0.2500 0.1875 -0.3125 -0.4375 -0.0625 -1.5625 0.1250 -1.3750 -0.8750 41 42 43 44 45 46 47 48 49 50 -0.8125 -1.1875 -1.3125 -2.5000 -1.0625 1.9375 0.7500 2.1875 1.6875 1.5625 51 52 53 54 55 56 57 58 59 60 3.9375 4.4375 4.1250 2.6250 1.1250 2.1875 2.8125 2.6875 1.5000 -1.0625 61 62 63 64 65 66 67 68 69 70 0.9375 0.7500 0.1875 3.6875 2.5625 -1.0625 0.4375 -0.8750 1.6250 1.1250 71 72 73 74 75 76 77 78 79 80 1.1875 0.8125 0.6875 3.5000 1.9375 -0.0625 0.7500 0.1875 -0.3125 -2.4375 81 82 83 84 85 86 87 88 89 90 -0.0625 0.4375 0.1250 -1.3750 -2.8750 -1.8125 -2.1875 -2.3125 -3.5000 -2.0625 91 92 93 94 95 96 97 98 99 100 0.9375 0.7500 0.1875 -0.3125 1.5625 0.9375 0.4375 0.1250 2.6250 2.1250 101 102 103 104 105 106 107 108 109 110 -1.8125 -2.1875 -2.3125 -0.5000 -0.0625 0.9375 1.7500 2.1875 1.6875 1.5625 111 112 113 114 115 116 117 118 119 120 1.9375 1.4375 1.1250 -0.3750 2.1250 0.1875 0.8125 1.6875 0.5000 -0.0625 121 122 123 124 125 126 127 128 129 130 -2.0625 -1.2500 -1.8125 -2.3125 -4.4375 -3.0625 -3.5625 -3.8750 -1.3750 -1.8750 131 132 133 134 135 136 137 138 139 140 0.1875 -0.1875 -0.3125 -1.5000 -0.0625 0.9375 -0.2500 1.1875 0.6875 0.5625 141 142 143 144 145 146 147 148 149 150 -1.0625 -0.5625 0.1250 -1.3750 -0.8750 1.1875 1.8125 -0.3125 2.5000 2.9375 151 152 153 154 155 156 157 158 159 160 -2.0625 -2.2500 -2.8125 -3.3125 -1.4375 -3.0625 -3.5625 -3.8750 -1.3750 -1.8750 161 162 163 164 165 166 167 168 169 170 1.1875 -0.1875 1.6875 0.5000 -0.0625 -2.0625 -3.2500 -1.8125 -2.3125 -2.4375 171 172 173 174 175 176 177 178 179 180 0.9375 1.4375 1.1250 -0.3750 -1.8750 -0.8125 -0.1875 0.6875 -0.5000 -1.0625 181 182 183 184 185 186 187 188 189 190 2.9375 3.7500 3.1875 2.6875 0.5625 3.9375 3.4375 3.1250 1.6250 4.1250 191 192 193 194 195 196 197 198 199 200 3.1875 3.8125 1.6875 4.5000 4.9375 1.9375 2.7500 3.1875 2.6875 2.5625 201 202 203 204 205 206 207 208 209 210 1.9375 1.4375 1.1250 3.6250 3.1250 1.1875 0.8125 0.6875 -0.5000 0.9375 211 212 213 214 215 216 217 218 219 220 -0.0625 -0.2500 -0.8125 -1.3125 0.5625 -3.0625 -1.5625 1.1250 -0.3750 0.1250 221 222 223 224 225 226 227 228 229 230 -1.8125 -1.1875 -0.3125 -1.5000 -2.0625 -1.0625 -2.2500 -0.8125 -1.3125 -1.4375 231 232 233 234 235 236 237 238 239 240 -1.0625 -0.5625 -0.8750 -2.3750 -3.8750 -0.8125 -1.1875 -1.3125 1.5000 -0.0625 > residuals(mixed) Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data. residuals(..., append = TRUE) will return data and residuals. 1 2 3 4 5 -2.333333e+00 -1.666667e+00 -3.666667e+00 -1.000000e+00 3.333333e-01 6 7 8 9 10 6.661338e-16 -1.000000e+00 2.775558e-17 -1.333333e+00 -1.000000e+00 11 12 13 14 15 -2.333333e+00 -2.000000e+00 -6.666667e-01 -2.000000e+00 -3.000000e+00 16 17 18 19 20 1.000000e+00 1.000000e+00 1.110223e-16 8.326673e-17 1.333333e+00 21 22 23 24 25 -6.666667e-01 -2.775558e-16 -2.000000e+00 6.666667e-01 2.000000e+00 26 27 28 29 30 -3.333333e-01 -6.666667e-01 -6.666667e-01 -6.666667e-01 6.666667e-01 31 32 33 34 35 -2.000000e+00 -1.000000e+00 -1.000000e+00 -1.000000e+00 -1.666667e+00 36 37 38 39 40 -1.666667e+00 -3.000000e+00 -1.000000e+00 -2.333333e+00 -2.000000e+00 41 42 43 44 45 -1.333333e+00 -1.666667e+00 -1.666667e+00 -1.666667e+00 -3.333333e-01 46 47 48 49 50 1.000000e+00 -2.775558e-16 1.000000e+00 1.000000e+00 3.333333e-01 51 52 53 54 55 2.333333e+00 3.000000e+00 3.000000e+00 1.666667e+00 8.326673e-17 56 57 58 59 60 1.666667e+00 2.333333e+00 2.333333e+00 2.333333e+00 -3.333333e-01 61 62 63 64 65 2.500000e-01 -2.500000e-01 -5.000000e-01 2.500000e+00 1.750000e+00 66 67 68 69 70 -1.500000e+00 -2.500000e-01 -1.000000e+00 1.000000e+00 5.000000e-01 71 72 73 74 75 1.750000e+00 1.500000e+00 1.250000e+00 3.500000e+00 2.000000e+00 76 77 78 79 80 -7.500000e-01 -2.500000e-01 -5.000000e-01 -1.500000e+00 -3.250000e+00 81 82 83 84 85 -5.000000e-01 -2.500000e-01 -2.893311e-16 -2.000000e+00 -3.500000e+00 86 87 88 89 90 -1.250000e+00 -1.500000e+00 -1.750000e+00 -3.500000e+00 -2.000000e+00 91 92 93 94 95 2.500000e-01 -2.500000e-01 -5.000000e-01 -1.500000e+00 7.500000e-01 96 97 98 99 100 5.000000e-01 -2.500000e-01 -2.893311e-16 2.000000e+00 1.500000e+00 101 102 103 104 105 -1.250000e+00 -1.500000e+00 -1.750000e+00 -5.000000e-01 5.551115e-16 106 107 108 109 110 2.500000e-01 7.500000e-01 1.500000e+00 5.000000e-01 7.500000e-01 111 112 113 114 115 1.500000e+00 7.500000e-01 1.000000e+00 -1.000000e+00 1.500000e+00 116 117 118 119 120 7.500000e-01 1.500000e+00 2.250000e+00 5.000000e-01 5.551115e-16 121 122 123 124 125 -3.333333e-01 3.333333e-01 3.333333e-01 -1.000000e+00 -2.666667e+00 126 127 128 129 130 -1.000000e+00 -1.000000e+00 -2.000000e+00 6.666667e-01 1.582068e-15 131 132 133 134 135 6.666667e-01 4.440892e-16 3.333333e-01 -1.000000e+00 1.110223e-16 136 137 138 139 140 2.666667e+00 1.333333e+00 3.333333e+00 2.000000e+00 2.333333e+00 141 142 143 144 145 1.000000e+00 2.000000e+00 2.000000e+00 6.666667e-01 1.000000e+00 146 147 148 149 150 1.666667e+00 2.000000e+00 3.333333e-01 3.000000e+00 3.000000e+00 151 152 153 154 155 4.440892e-16 5.000000e-01 -5.000000e-01 -5.000000e-01 5.000000e-01 156 157 158 159 160 -2.000000e+00 -2.500000e+00 -2.500000e+00 -5.000000e-01 -3.330669e-16 161 162 163 164 165 1.000000e+00 -4.440892e-16 5.000000e-01 5.000000e-01 5.000000e-01 166 167 168 169 170 -1.110223e-16 -5.000000e-01 5.000000e-01 5.000000e-01 -5.000000e-01 171 172 173 174 175 2.000000e+00 2.500000e+00 2.500000e+00 5.000000e-01 4.718448e-16 176 177 178 179 180 -1.000000e+00 9.436896e-16 -5.000000e-01 -5.000000e-01 -5.000000e-01 181 182 183 184 185 5.000000e-01 5.000000e-01 -1.110223e-15 -8.049117e-16 -1.000000e+00 186 187 188 189 190 1.000000e+00 1.000000e+00 1.000000e+00 -1.000000e+00 5.000000e-01 191 192 193 194 195 1.000000e+00 1.500000e+00 5.000000e-01 2.500000e+00 2.000000e+00 196 197 198 199 200 -5.000000e-01 -5.000000e-01 2.220446e-16 -5.828671e-16 1.000000e+00 201 202 203 204 205 -1.000000e+00 -1.000000e+00 -1.000000e+00 1.000000e+00 -5.000000e-01 206 207 208 209 210 -1.000000e+00 -1.500000e+00 -5.000000e-01 -2.500000e+00 -2.000000e+00 211 212 213 214 215 5.000000e-01 1.000000e+00 -2.220446e-16 -8.049117e-16 1.000000e+00 216 217 218 219 220 -1.000000e+00 -5.000000e-01 1.000000e+00 1.000000e+00 2.000000e+00 221 222 223 224 225 -5.000000e-01 -6.106227e-16 5.000000e-01 -1.500000e+00 -1.000000e+00 226 227 228 229 230 -5.000000e-01 -1.000000e+00 -2.220446e-16 -8.049117e-16 -1.000000e+00 231 232 233 234 235 1.000000e+00 5.000000e-01 -1.000000e+00 -1.000000e+00 -2.000000e+00 236 237 238 239 240 5.000000e-01 -6.106227e-16 -5.000000e-01 1.500000e+00 1.000000e+00 > residuals(between) Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data. residuals(..., append = TRUE) will return data and residuals. 1 2 3 4 5 6 7 -1.4444444 -0.4444444 1.8888889 -0.3333333 0.3333333 0.6666667 -0.6666667 8 9 10 11 12 13 14 0.1666667 -0.1666667 0.1111111 -1.5555556 1.4444444 0.8333333 -1.5000000 15 16 -0.1666667 0.8333333 > > ## Get residuals plus data used for fitting: > residuals(within, append = TRUE) id phase hour value .residuals 1 1 fup X1 2 -4.0625 2 1 fup X2 3 -3.2500 3 1 fup X3 2 -5.8125 4 1 fup X4 4 -2.3125 5 1 fup X5 4 -1.4375 6 1 post X1 3 -2.0625 7 1 post X2 2 -3.5625 8 1 post X3 5 -1.8750 9 1 post X4 3 -3.3750 10 1 post X5 2 -2.8750 11 1 pre X1 1 -2.8125 12 1 pre X2 2 -2.1875 13 1 pre X3 4 -1.3125 14 1 pre X4 2 -2.5000 15 1 pre X5 1 -3.0625 16 10 fup X1 8 1.9375 17 10 fup X2 8 1.7500 18 10 fup X3 9 1.1875 19 10 fup X4 7 0.6875 20 10 fup X5 8 2.5625 21 10 post X1 6 0.9375 22 10 post X2 7 1.4375 23 10 post X3 6 -0.8750 24 10 post X4 8 1.6250 25 10 post X5 8 3.1250 26 10 pre X1 4 0.1875 27 10 pre X2 4 -0.1875 28 10 pre X3 5 -0.3125 29 10 pre X4 3 -1.5000 30 10 pre X5 4 -0.0625 31 11 fup X1 5 -1.0625 32 11 fup X2 6 -0.2500 33 11 fup X3 8 0.1875 34 11 fup X4 6 -0.3125 35 11 fup X5 5 -0.4375 36 11 post X1 5 -0.0625 37 11 post X2 4 -1.5625 38 11 post X3 7 0.1250 39 11 post X4 5 -1.3750 40 11 post X5 4 -0.8750 41 11 pre X1 3 -0.8125 42 11 pre X2 3 -1.1875 43 11 pre X3 4 -1.3125 44 11 pre X4 2 -2.5000 45 11 pre X5 3 -1.0625 46 12 fup X1 8 1.9375 47 12 fup X2 7 0.7500 48 12 fup X3 10 2.1875 49 12 fup X4 8 1.6875 50 12 fup X5 7 1.5625 51 12 post X1 9 3.9375 52 12 post X2 10 4.4375 53 12 post X3 11 4.1250 54 12 post X4 9 2.6250 55 12 post X5 6 1.1250 56 12 pre X1 6 2.1875 57 12 pre X2 7 2.8125 58 12 pre X3 8 2.6875 59 12 pre X4 6 1.5000 60 12 pre X5 3 -1.0625 61 13 fup X1 7 0.9375 62 13 fup X2 7 0.7500 63 13 fup X3 8 0.1875 64 13 fup X4 10 3.6875 65 13 fup X5 8 2.5625 66 13 post X1 4 -1.0625 67 13 post X2 6 0.4375 68 13 post X3 6 -0.8750 69 13 post X4 8 1.6250 70 13 post X5 6 1.1250 71 13 pre X1 5 1.1875 72 13 pre X2 5 0.8125 73 13 pre X3 6 0.6875 74 13 pre X4 8 3.5000 75 13 pre X5 6 1.9375 76 14 fup X1 6 -0.0625 77 14 fup X2 7 0.7500 78 14 fup X3 8 0.1875 79 14 fup X4 6 -0.3125 80 14 fup X5 3 -2.4375 81 14 post X1 5 -0.0625 82 14 post X2 6 0.4375 83 14 post X3 7 0.1250 84 14 post X4 5 -1.3750 85 14 post X5 2 -2.8750 86 14 pre X1 2 -1.8125 87 14 pre X2 2 -2.1875 88 14 pre X3 3 -2.3125 89 14 pre X4 1 -3.5000 90 14 pre X5 2 -2.0625 91 15 fup X1 7 0.9375 92 15 fup X2 7 0.7500 93 15 fup X3 8 0.1875 94 15 fup X4 6 -0.3125 95 15 fup X5 7 1.5625 96 15 post X1 6 0.9375 97 15 post X2 6 0.4375 98 15 post X3 7 0.1250 99 15 post X4 9 2.6250 100 15 post X5 7 2.1250 101 15 pre X1 2 -1.8125 102 15 pre X2 2 -2.1875 103 15 pre X3 3 -2.3125 104 15 pre X4 4 -0.5000 105 15 pre X5 4 -0.0625 106 16 fup X1 7 0.9375 107 16 fup X2 8 1.7500 108 16 fup X3 10 2.1875 109 16 fup X4 8 1.6875 110 16 fup X5 7 1.5625 111 16 post X1 7 1.9375 112 16 post X2 7 1.4375 113 16 post X3 8 1.1250 114 16 post X4 6 -0.3750 115 16 post X5 7 2.1250 116 16 pre X1 4 0.1875 117 16 pre X2 5 0.8125 118 16 pre X3 7 1.6875 119 16 pre X4 5 0.5000 120 16 pre X5 4 -0.0625 121 2 fup X1 4 -2.0625 122 2 fup X2 5 -1.2500 123 2 fup X3 6 -1.8125 124 2 fup X4 4 -2.3125 125 2 fup X5 1 -4.4375 126 2 post X1 2 -3.0625 127 2 post X2 2 -3.5625 128 2 post X3 3 -3.8750 129 2 post X4 5 -1.3750 130 2 post X5 3 -1.8750 131 2 pre X1 4 0.1875 132 2 pre X2 4 -0.1875 133 2 pre X3 5 -0.3125 134 2 pre X4 3 -1.5000 135 2 pre X5 4 -0.0625 136 3 fup X1 7 0.9375 137 3 fup X2 6 -0.2500 138 3 fup X3 9 1.1875 139 3 fup X4 7 0.6875 140 3 fup X5 6 0.5625 141 3 post X1 4 -1.0625 142 3 post X2 5 -0.5625 143 3 post X3 7 0.1250 144 3 post X4 5 -1.3750 145 3 post X5 4 -0.8750 146 3 pre X1 5 1.1875 147 3 pre X2 6 1.8125 148 3 pre X3 5 -0.3125 149 3 pre X4 7 2.5000 150 3 pre X5 7 2.9375 151 4 fup X1 4 -2.0625 152 4 fup X2 4 -2.2500 153 4 fup X3 5 -2.8125 154 4 fup X4 3 -3.3125 155 4 fup X5 4 -1.4375 156 4 post X1 2 -3.0625 157 4 post X2 2 -3.5625 158 4 post X3 3 -3.8750 159 4 post X4 5 -1.3750 160 4 post X5 3 -1.8750 161 4 pre X1 5 1.1875 162 4 pre X2 4 -0.1875 163 4 pre X3 7 1.6875 164 4 pre X4 5 0.5000 165 4 pre X5 4 -0.0625 166 5 fup X1 4 -2.0625 167 5 fup X2 3 -3.2500 168 5 fup X3 6 -1.8125 169 5 fup X4 4 -2.3125 170 5 fup X5 3 -2.4375 171 5 post X1 6 0.9375 172 5 post X2 7 1.4375 173 5 post X3 8 1.1250 174 5 post X4 6 -0.3750 175 5 post X5 3 -1.8750 176 5 pre X1 3 -0.8125 177 5 pre X2 4 -0.1875 178 5 pre X3 6 0.6875 179 5 pre X4 4 -0.5000 180 5 pre X5 3 -1.0625 181 6 fup X1 9 2.9375 182 6 fup X2 10 3.7500 183 6 fup X3 11 3.1875 184 6 fup X4 9 2.6875 185 6 fup X5 6 0.5625 186 6 post X1 9 3.9375 187 6 post X2 9 3.4375 188 6 post X3 10 3.1250 189 6 post X4 8 1.6250 190 6 post X5 9 4.1250 191 6 pre X1 7 3.1875 192 6 pre X2 8 3.8125 193 6 pre X3 7 1.6875 194 6 pre X4 9 4.5000 195 6 pre X5 9 4.9375 196 7 fup X1 8 1.9375 197 7 fup X2 9 2.7500 198 7 fup X3 11 3.1875 199 7 fup X4 9 2.6875 200 7 fup X5 8 2.5625 201 7 post X1 7 1.9375 202 7 post X2 7 1.4375 203 7 post X3 8 1.1250 204 7 post X4 10 3.6250 205 7 post X5 8 3.1250 206 7 pre X1 5 1.1875 207 7 pre X2 5 0.8125 208 7 pre X3 6 0.6875 209 7 pre X4 4 -0.5000 210 7 pre X5 5 0.9375 211 8 fup X1 6 -0.0625 212 8 fup X2 6 -0.2500 213 8 fup X3 7 -0.8125 214 8 fup X4 5 -1.3125 215 8 fup X5 6 0.5625 216 8 post X1 2 -3.0625 217 8 post X2 4 -1.5625 218 8 post X3 8 1.1250 219 8 post X4 6 -0.3750 220 8 post X5 5 0.1250 221 8 pre X1 2 -1.8125 222 8 pre X2 3 -1.1875 223 8 pre X3 5 -0.3125 224 8 pre X4 3 -1.5000 225 8 pre X5 2 -2.0625 226 9 fup X1 5 -1.0625 227 9 fup X2 4 -2.2500 228 9 fup X3 7 -0.8125 229 9 fup X4 5 -1.3125 230 9 fup X5 4 -1.4375 231 9 post X1 4 -1.0625 232 9 post X2 5 -0.5625 233 9 post X3 6 -0.8750 234 9 post X4 4 -2.3750 235 9 post X5 1 -3.8750 236 9 pre X1 3 -0.8125 237 9 pre X2 3 -1.1875 238 9 pre X3 4 -1.3125 239 9 pre X4 6 1.5000 240 9 pre X5 4 -0.0625 > residuals(mixed, append = TRUE) id phase hour treatment gender value .residuals 1 1 fup X1 control M 2 -2.333333e+00 2 1 fup X2 control M 3 -1.666667e+00 3 1 fup X3 control M 2 -3.666667e+00 4 1 fup X4 control M 4 -1.000000e+00 5 1 fup X5 control M 4 3.333333e-01 6 1 post X1 control M 3 6.661338e-16 7 1 post X2 control M 2 -1.000000e+00 8 1 post X3 control M 5 2.775558e-17 9 1 post X4 control M 3 -1.333333e+00 10 1 post X5 control M 2 -1.000000e+00 11 1 pre X1 control M 1 -2.333333e+00 12 1 pre X2 control M 2 -2.000000e+00 13 1 pre X3 control M 4 -6.666667e-01 14 1 pre X4 control M 2 -2.000000e+00 15 1 pre X5 control M 1 -3.000000e+00 16 10 fup X1 B M 8 1.000000e+00 17 10 fup X2 B M 8 1.000000e+00 18 10 fup X3 B M 9 1.110223e-16 19 10 fup X4 B M 7 8.326673e-17 20 10 fup X5 B M 8 1.333333e+00 21 10 post X1 B M 6 -6.666667e-01 22 10 post X2 B M 7 -2.775558e-16 23 10 post X3 B M 6 -2.000000e+00 24 10 post X4 B M 8 6.666667e-01 25 10 post X5 B M 8 2.000000e+00 26 10 pre X1 B M 4 -3.333333e-01 27 10 pre X2 B M 4 -6.666667e-01 28 10 pre X3 B M 5 -6.666667e-01 29 10 pre X4 B M 3 -6.666667e-01 30 10 pre X5 B M 4 6.666667e-01 31 11 fup X1 B M 5 -2.000000e+00 32 11 fup X2 B M 6 -1.000000e+00 33 11 fup X3 B M 8 -1.000000e+00 34 11 fup X4 B M 6 -1.000000e+00 35 11 fup X5 B M 5 -1.666667e+00 36 11 post X1 B M 5 -1.666667e+00 37 11 post X2 B M 4 -3.000000e+00 38 11 post X3 B M 7 -1.000000e+00 39 11 post X4 B M 5 -2.333333e+00 40 11 post X5 B M 4 -2.000000e+00 41 11 pre X1 B M 3 -1.333333e+00 42 11 pre X2 B M 3 -1.666667e+00 43 11 pre X3 B M 4 -1.666667e+00 44 11 pre X4 B M 2 -1.666667e+00 45 11 pre X5 B M 3 -3.333333e-01 46 12 fup X1 B M 8 1.000000e+00 47 12 fup X2 B M 7 -2.775558e-16 48 12 fup X3 B M 10 1.000000e+00 49 12 fup X4 B M 8 1.000000e+00 50 12 fup X5 B M 7 3.333333e-01 51 12 post X1 B M 9 2.333333e+00 52 12 post X2 B M 10 3.000000e+00 53 12 post X3 B M 11 3.000000e+00 54 12 post X4 B M 9 1.666667e+00 55 12 post X5 B M 6 8.326673e-17 56 12 pre X1 B M 6 1.666667e+00 57 12 pre X2 B M 7 2.333333e+00 58 12 pre X3 B M 8 2.333333e+00 59 12 pre X4 B M 6 2.333333e+00 60 12 pre X5 B M 3 -3.333333e-01 61 13 fup X1 B F 7 2.500000e-01 62 13 fup X2 B F 7 -2.500000e-01 63 13 fup X3 B F 8 -5.000000e-01 64 13 fup X4 B F 10 2.500000e+00 65 13 fup X5 B F 8 1.750000e+00 66 13 post X1 B F 4 -1.500000e+00 67 13 post X2 B F 6 -2.500000e-01 68 13 post X3 B F 6 -1.000000e+00 69 13 post X4 B F 8 1.000000e+00 70 13 post X5 B F 6 5.000000e-01 71 13 pre X1 B F 5 1.750000e+00 72 13 pre X2 B F 5 1.500000e+00 73 13 pre X3 B F 6 1.250000e+00 74 13 pre X4 B F 8 3.500000e+00 75 13 pre X5 B F 6 2.000000e+00 76 14 fup X1 B F 6 -7.500000e-01 77 14 fup X2 B F 7 -2.500000e-01 78 14 fup X3 B F 8 -5.000000e-01 79 14 fup X4 B F 6 -1.500000e+00 80 14 fup X5 B F 3 -3.250000e+00 81 14 post X1 B F 5 -5.000000e-01 82 14 post X2 B F 6 -2.500000e-01 83 14 post X3 B F 7 -2.893311e-16 84 14 post X4 B F 5 -2.000000e+00 85 14 post X5 B F 2 -3.500000e+00 86 14 pre X1 B F 2 -1.250000e+00 87 14 pre X2 B F 2 -1.500000e+00 88 14 pre X3 B F 3 -1.750000e+00 89 14 pre X4 B F 1 -3.500000e+00 90 14 pre X5 B F 2 -2.000000e+00 91 15 fup X1 B F 7 2.500000e-01 92 15 fup X2 B F 7 -2.500000e-01 93 15 fup X3 B F 8 -5.000000e-01 94 15 fup X4 B F 6 -1.500000e+00 95 15 fup X5 B F 7 7.500000e-01 96 15 post X1 B F 6 5.000000e-01 97 15 post X2 B F 6 -2.500000e-01 98 15 post X3 B F 7 -2.893311e-16 99 15 post X4 B F 9 2.000000e+00 100 15 post X5 B F 7 1.500000e+00 101 15 pre X1 B F 2 -1.250000e+00 102 15 pre X2 B F 2 -1.500000e+00 103 15 pre X3 B F 3 -1.750000e+00 104 15 pre X4 B F 4 -5.000000e-01 105 15 pre X5 B F 4 5.551115e-16 106 16 fup X1 B F 7 2.500000e-01 107 16 fup X2 B F 8 7.500000e-01 108 16 fup X3 B F 10 1.500000e+00 109 16 fup X4 B F 8 5.000000e-01 110 16 fup X5 B F 7 7.500000e-01 111 16 post X1 B F 7 1.500000e+00 112 16 post X2 B F 7 7.500000e-01 113 16 post X3 B F 8 1.000000e+00 114 16 post X4 B F 6 -1.000000e+00 115 16 post X5 B F 7 1.500000e+00 116 16 pre X1 B F 4 7.500000e-01 117 16 pre X2 B F 5 1.500000e+00 118 16 pre X3 B F 7 2.250000e+00 119 16 pre X4 B F 5 5.000000e-01 120 16 pre X5 B F 4 5.551115e-16 121 2 fup X1 control M 4 -3.333333e-01 122 2 fup X2 control M 5 3.333333e-01 123 2 fup X3 control M 6 3.333333e-01 124 2 fup X4 control M 4 -1.000000e+00 125 2 fup X5 control M 1 -2.666667e+00 126 2 post X1 control M 2 -1.000000e+00 127 2 post X2 control M 2 -1.000000e+00 128 2 post X3 control M 3 -2.000000e+00 129 2 post X4 control M 5 6.666667e-01 130 2 post X5 control M 3 1.582068e-15 131 2 pre X1 control M 4 6.666667e-01 132 2 pre X2 control M 4 4.440892e-16 133 2 pre X3 control M 5 3.333333e-01 134 2 pre X4 control M 3 -1.000000e+00 135 2 pre X5 control M 4 1.110223e-16 136 3 fup X1 control M 7 2.666667e+00 137 3 fup X2 control M 6 1.333333e+00 138 3 fup X3 control M 9 3.333333e+00 139 3 fup X4 control M 7 2.000000e+00 140 3 fup X5 control M 6 2.333333e+00 141 3 post X1 control M 4 1.000000e+00 142 3 post X2 control M 5 2.000000e+00 143 3 post X3 control M 7 2.000000e+00 144 3 post X4 control M 5 6.666667e-01 145 3 post X5 control M 4 1.000000e+00 146 3 pre X1 control M 5 1.666667e+00 147 3 pre X2 control M 6 2.000000e+00 148 3 pre X3 control M 5 3.333333e-01 149 3 pre X4 control M 7 3.000000e+00 150 3 pre X5 control M 7 3.000000e+00 151 4 fup X1 control F 4 4.440892e-16 152 4 fup X2 control F 4 5.000000e-01 153 4 fup X3 control F 5 -5.000000e-01 154 4 fup X4 control F 3 -5.000000e-01 155 4 fup X5 control F 4 5.000000e-01 156 4 post X1 control F 2 -2.000000e+00 157 4 post X2 control F 2 -2.500000e+00 158 4 post X3 control F 3 -2.500000e+00 159 4 post X4 control F 5 -5.000000e-01 160 4 post X5 control F 3 -3.330669e-16 161 4 pre X1 control F 5 1.000000e+00 162 4 pre X2 control F 4 -4.440892e-16 163 4 pre X3 control F 7 5.000000e-01 164 4 pre X4 control F 5 5.000000e-01 165 4 pre X5 control F 4 5.000000e-01 166 5 fup X1 control F 4 -1.110223e-16 167 5 fup X2 control F 3 -5.000000e-01 168 5 fup X3 control F 6 5.000000e-01 169 5 fup X4 control F 4 5.000000e-01 170 5 fup X5 control F 3 -5.000000e-01 171 5 post X1 control F 6 2.000000e+00 172 5 post X2 control F 7 2.500000e+00 173 5 post X3 control F 8 2.500000e+00 174 5 post X4 control F 6 5.000000e-01 175 5 post X5 control F 3 4.718448e-16 176 5 pre X1 control F 3 -1.000000e+00 177 5 pre X2 control F 4 9.436896e-16 178 5 pre X3 control F 6 -5.000000e-01 179 5 pre X4 control F 4 -5.000000e-01 180 5 pre X5 control F 3 -5.000000e-01 181 6 fup X1 A M 9 5.000000e-01 182 6 fup X2 A M 10 5.000000e-01 183 6 fup X3 A M 11 -1.110223e-15 184 6 fup X4 A M 9 -8.049117e-16 185 6 fup X5 A M 6 -1.000000e+00 186 6 post X1 A M 9 1.000000e+00 187 6 post X2 A M 9 1.000000e+00 188 6 post X3 A M 10 1.000000e+00 189 6 post X4 A M 8 -1.000000e+00 190 6 post X5 A M 9 5.000000e-01 191 6 pre X1 A M 7 1.000000e+00 192 6 pre X2 A M 8 1.500000e+00 193 6 pre X3 A M 7 5.000000e-01 194 6 pre X4 A M 9 2.500000e+00 195 6 pre X5 A M 9 2.000000e+00 196 7 fup X1 A M 8 -5.000000e-01 197 7 fup X2 A M 9 -5.000000e-01 198 7 fup X3 A M 11 2.220446e-16 199 7 fup X4 A M 9 -5.828671e-16 200 7 fup X5 A M 8 1.000000e+00 201 7 post X1 A M 7 -1.000000e+00 202 7 post X2 A M 7 -1.000000e+00 203 7 post X3 A M 8 -1.000000e+00 204 7 post X4 A M 10 1.000000e+00 205 7 post X5 A M 8 -5.000000e-01 206 7 pre X1 A M 5 -1.000000e+00 207 7 pre X2 A M 5 -1.500000e+00 208 7 pre X3 A M 6 -5.000000e-01 209 7 pre X4 A M 4 -2.500000e+00 210 7 pre X5 A M 5 -2.000000e+00 211 8 fup X1 A F 6 5.000000e-01 212 8 fup X2 A F 6 1.000000e+00 213 8 fup X3 A F 7 -2.220446e-16 214 8 fup X4 A F 5 -8.049117e-16 215 8 fup X5 A F 6 1.000000e+00 216 8 post X1 A F 2 -1.000000e+00 217 8 post X2 A F 4 -5.000000e-01 218 8 post X3 A F 8 1.000000e+00 219 8 post X4 A F 6 1.000000e+00 220 8 post X5 A F 5 2.000000e+00 221 8 pre X1 A F 2 -5.000000e-01 222 8 pre X2 A F 3 -6.106227e-16 223 8 pre X3 A F 5 5.000000e-01 224 8 pre X4 A F 3 -1.500000e+00 225 8 pre X5 A F 2 -1.000000e+00 226 9 fup X1 A F 5 -5.000000e-01 227 9 fup X2 A F 4 -1.000000e+00 228 9 fup X3 A F 7 -2.220446e-16 229 9 fup X4 A F 5 -8.049117e-16 230 9 fup X5 A F 4 -1.000000e+00 231 9 post X1 A F 4 1.000000e+00 232 9 post X2 A F 5 5.000000e-01 233 9 post X3 A F 6 -1.000000e+00 234 9 post X4 A F 4 -1.000000e+00 235 9 post X5 A F 1 -2.000000e+00 236 9 pre X1 A F 3 5.000000e-01 237 9 pre X2 A F 3 -6.106227e-16 238 9 pre X3 A F 4 -5.000000e-01 239 9 pre X4 A F 6 1.500000e+00 240 9 pre X5 A F 4 1.000000e+00 > residuals(between, append = TRUE) id treatment gender value .residuals 1 1 control M 2.666667 -1.4444444 2 2 control M 3.666667 -0.4444444 3 3 control M 6.000000 1.8888889 4 4 control F 4.000000 -0.3333333 5 5 control F 4.666667 0.3333333 6 6 A M 8.666667 0.6666667 7 7 A M 7.333333 -0.6666667 8 8 A F 4.666667 0.1666667 9 9 A F 4.333333 -0.1666667 10 10 B M 6.333333 0.1111111 11 11 B M 4.666667 -1.5555556 12 12 B M 7.666667 1.4444444 13 13 B F 6.666667 0.8333333 14 14 B F 4.333333 -1.5000000 15 15 B F 5.666667 -0.1666667 16 16 B F 6.666667 0.8333333 > > ### in case data is correctly ordered before fitting, this message is not shown > > ## between data: > obk2 <- aggregate(value ~ gender + treatment + id , data = obk.long, FUN = mean) > between2 <- aov_car(value ~ treatment*gender + Error(id), data = obk2) Contrasts set to contr.sum for the following variables: treatment, gender > > residuals(between2) ## no message 1 2 3 4 5 6 7 -1.4444444 -0.4444444 1.8888889 -0.3333333 0.3333333 0.6666667 -0.6666667 8 9 10 11 12 13 14 0.1666667 -0.1666667 0.1111111 -1.5555556 1.4444444 0.8333333 -1.5000000 15 16 -0.1666667 0.8333333 > all.equal(obk2, between2$data$long[,colnames(obk2)]) ## TRUE [1] TRUE > > # Therefore okay: > obk2$residuals <- residuals(between2) > > ## within data > obk3 <- obk.long[with(obk.long, order(id, phase, hour)), ] > within2 <- aov_car(value ~ 1 + Error(id/(phase*hour)), data = obk3) > residuals(within2) ## no message, because order is correct 1 2 3 4 5 6 7 8 9 10 -4.0625 -3.2500 -5.8125 -2.3125 -1.4375 -2.0625 -3.5625 -1.8750 -3.3750 -2.8750 11 12 13 14 15 16 17 18 19 20 -2.8125 -2.1875 -1.3125 -2.5000 -3.0625 1.9375 1.7500 1.1875 0.6875 2.5625 21 22 23 24 25 26 27 28 29 30 0.9375 1.4375 -0.8750 1.6250 3.1250 0.1875 -0.1875 -0.3125 -1.5000 -0.0625 31 32 33 34 35 36 37 38 39 40 -1.0625 -0.2500 0.1875 -0.3125 -0.4375 -0.0625 -1.5625 0.1250 -1.3750 -0.8750 41 42 43 44 45 46 47 48 49 50 -0.8125 -1.1875 -1.3125 -2.5000 -1.0625 1.9375 0.7500 2.1875 1.6875 1.5625 51 52 53 54 55 56 57 58 59 60 3.9375 4.4375 4.1250 2.6250 1.1250 2.1875 2.8125 2.6875 1.5000 -1.0625 61 62 63 64 65 66 67 68 69 70 0.9375 0.7500 0.1875 3.6875 2.5625 -1.0625 0.4375 -0.8750 1.6250 1.1250 71 72 73 74 75 76 77 78 79 80 1.1875 0.8125 0.6875 3.5000 1.9375 -0.0625 0.7500 0.1875 -0.3125 -2.4375 81 82 83 84 85 86 87 88 89 90 -0.0625 0.4375 0.1250 -1.3750 -2.8750 -1.8125 -2.1875 -2.3125 -3.5000 -2.0625 91 92 93 94 95 96 97 98 99 100 0.9375 0.7500 0.1875 -0.3125 1.5625 0.9375 0.4375 0.1250 2.6250 2.1250 101 102 103 104 105 106 107 108 109 110 -1.8125 -2.1875 -2.3125 -0.5000 -0.0625 0.9375 1.7500 2.1875 1.6875 1.5625 111 112 113 114 115 116 117 118 119 120 1.9375 1.4375 1.1250 -0.3750 2.1250 0.1875 0.8125 1.6875 0.5000 -0.0625 121 122 123 124 125 126 127 128 129 130 -2.0625 -1.2500 -1.8125 -2.3125 -4.4375 -3.0625 -3.5625 -3.8750 -1.3750 -1.8750 131 132 133 134 135 136 137 138 139 140 0.1875 -0.1875 -0.3125 -1.5000 -0.0625 0.9375 -0.2500 1.1875 0.6875 0.5625 141 142 143 144 145 146 147 148 149 150 -1.0625 -0.5625 0.1250 -1.3750 -0.8750 1.1875 1.8125 -0.3125 2.5000 2.9375 151 152 153 154 155 156 157 158 159 160 -2.0625 -2.2500 -2.8125 -3.3125 -1.4375 -3.0625 -3.5625 -3.8750 -1.3750 -1.8750 161 162 163 164 165 166 167 168 169 170 1.1875 -0.1875 1.6875 0.5000 -0.0625 -2.0625 -3.2500 -1.8125 -2.3125 -2.4375 171 172 173 174 175 176 177 178 179 180 0.9375 1.4375 1.1250 -0.3750 -1.8750 -0.8125 -0.1875 0.6875 -0.5000 -1.0625 181 182 183 184 185 186 187 188 189 190 2.9375 3.7500 3.1875 2.6875 0.5625 3.9375 3.4375 3.1250 1.6250 4.1250 191 192 193 194 195 196 197 198 199 200 3.1875 3.8125 1.6875 4.5000 4.9375 1.9375 2.7500 3.1875 2.6875 2.5625 201 202 203 204 205 206 207 208 209 210 1.9375 1.4375 1.1250 3.6250 3.1250 1.1875 0.8125 0.6875 -0.5000 0.9375 211 212 213 214 215 216 217 218 219 220 -0.0625 -0.2500 -0.8125 -1.3125 0.5625 -3.0625 -1.5625 1.1250 -0.3750 0.1250 221 222 223 224 225 226 227 228 229 230 -1.8125 -1.1875 -0.3125 -1.5000 -2.0625 -1.0625 -2.2500 -0.8125 -1.3125 -1.4375 231 232 233 234 235 236 237 238 239 240 -1.0625 -0.5625 -0.8750 -2.3750 -3.8750 -0.8125 -1.1875 -1.3125 1.5000 -0.0625 > # Therefore okay: > obk3$residuals <- residuals(within2) > > ## Same for fitted values: > # (show message) > fitted(within) Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data. fitted(..., append = TRUE) will return data and fitted values. 1 2 3 4 5 6 7 8 9 10 11 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 12 13 14 15 16 17 18 19 20 21 22 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 23 24 25 26 27 28 29 30 31 32 33 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 34 35 36 37 38 39 40 41 42 43 44 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 45 46 47 48 49 50 51 52 53 54 55 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 56 57 58 59 60 61 62 63 64 65 66 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 67 68 69 70 71 72 73 74 75 76 77 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 78 79 80 81 82 83 84 85 86 87 88 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 89 90 91 92 93 94 95 96 97 98 99 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 100 101 102 103 104 105 106 107 108 109 110 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 111 112 113 114 115 116 117 118 119 120 121 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 122 123 124 125 126 127 128 129 130 131 132 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 133 134 135 136 137 138 139 140 141 142 143 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 144 145 146 147 148 149 150 151 152 153 154 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 155 156 157 158 159 160 161 162 163 164 165 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 166 167 168 169 170 171 172 173 174 175 176 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 177 178 179 180 181 182 183 184 185 186 187 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 188 189 190 191 192 193 194 195 196 197 198 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 199 200 201 202 203 204 205 206 207 208 209 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 210 211 212 213 214 215 216 217 218 219 220 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 221 222 223 224 225 226 227 228 229 230 231 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 232 233 234 235 236 237 238 239 240 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 > fitted(mixed) Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data. fitted(..., append = TRUE) will return data and fitted values. 1 2 3 4 5 6 7 8 4.333333 4.666667 5.666667 5.000000 3.666667 3.000000 3.000000 5.000000 9 10 11 12 13 14 15 16 4.333333 3.000000 3.333333 4.000000 4.666667 4.000000 4.000000 7.000000 17 18 19 20 21 22 23 24 7.000000 9.000000 7.000000 6.666667 6.666667 7.000000 8.000000 7.333333 25 26 27 28 29 30 31 32 6.000000 4.333333 4.666667 5.666667 3.666667 3.333333 7.000000 7.000000 33 34 35 36 37 38 39 40 9.000000 7.000000 6.666667 6.666667 7.000000 8.000000 7.333333 6.000000 41 42 43 44 45 46 47 48 4.333333 4.666667 5.666667 3.666667 3.333333 7.000000 7.000000 9.000000 49 50 51 52 53 54 55 56 7.000000 6.666667 6.666667 7.000000 8.000000 7.333333 6.000000 4.333333 57 58 59 60 61 62 63 64 4.666667 5.666667 3.666667 3.333333 6.750000 7.250000 8.500000 7.500000 65 66 67 68 69 70 71 72 6.250000 5.500000 6.250000 7.000000 7.000000 5.500000 3.250000 3.500000 73 74 75 76 77 78 79 80 4.750000 4.500000 4.000000 6.750000 7.250000 8.500000 7.500000 6.250000 81 82 83 84 85 86 87 88 5.500000 6.250000 7.000000 7.000000 5.500000 3.250000 3.500000 4.750000 89 90 91 92 93 94 95 96 4.500000 4.000000 6.750000 7.250000 8.500000 7.500000 6.250000 5.500000 97 98 99 100 101 102 103 104 6.250000 7.000000 7.000000 5.500000 3.250000 3.500000 4.750000 4.500000 105 106 107 108 109 110 111 112 4.000000 6.750000 7.250000 8.500000 7.500000 6.250000 5.500000 6.250000 113 114 115 116 117 118 119 120 7.000000 7.000000 5.500000 3.250000 3.500000 4.750000 4.500000 4.000000 121 122 123 124 125 126 127 128 4.333333 4.666667 5.666667 5.000000 3.666667 3.000000 3.000000 5.000000 129 130 131 132 133 134 135 136 4.333333 3.000000 3.333333 4.000000 4.666667 4.000000 4.000000 4.333333 137 138 139 140 141 142 143 144 4.666667 5.666667 5.000000 3.666667 3.000000 3.000000 5.000000 4.333333 145 146 147 148 149 150 151 152 3.000000 3.333333 4.000000 4.666667 4.000000 4.000000 4.000000 3.500000 153 154 155 156 157 158 159 160 5.500000 3.500000 3.500000 4.000000 4.500000 5.500000 5.500000 3.000000 161 162 163 164 165 166 167 168 4.000000 4.000000 6.500000 4.500000 3.500000 4.000000 3.500000 5.500000 169 170 171 172 173 174 175 176 3.500000 3.500000 4.000000 4.500000 5.500000 5.500000 3.000000 4.000000 177 178 179 180 181 182 183 184 4.000000 6.500000 4.500000 3.500000 8.500000 9.500000 11.000000 9.000000 185 186 187 188 189 190 191 192 7.000000 8.000000 8.000000 9.000000 9.000000 8.500000 6.000000 6.500000 193 194 195 196 197 198 199 200 6.500000 6.500000 7.000000 8.500000 9.500000 11.000000 9.000000 7.000000 201 202 203 204 205 206 207 208 8.000000 8.000000 9.000000 9.000000 8.500000 6.000000 6.500000 6.500000 209 210 211 212 213 214 215 216 6.500000 7.000000 5.500000 5.000000 7.000000 5.000000 5.000000 3.000000 217 218 219 220 221 222 223 224 4.500000 7.000000 5.000000 3.000000 2.500000 3.000000 4.500000 4.500000 225 226 227 228 229 230 231 232 3.000000 5.500000 5.000000 7.000000 5.000000 5.000000 3.000000 4.500000 233 234 235 236 237 238 239 240 7.000000 5.000000 3.000000 2.500000 3.000000 4.500000 4.500000 3.000000 > fitted(between) Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data. fitted(..., append = TRUE) will return data and fitted values. 1 2 3 4 5 6 7 8 4.111111 4.111111 4.111111 4.333333 4.333333 8.000000 8.000000 4.500000 9 10 11 12 13 14 15 16 4.500000 6.222222 6.222222 6.222222 5.833333 5.833333 5.833333 5.833333 > > ## Get fitted values plus data used for fitting: > fitted(within, append = TRUE) id phase hour value .fitted 1 1 fup X1 2 6.0625 2 1 fup X2 3 6.2500 3 1 fup X3 2 7.8125 4 1 fup X4 4 6.3125 5 1 fup X5 4 5.4375 6 1 post X1 3 5.0625 7 1 post X2 2 5.5625 8 1 post X3 5 6.8750 9 1 post X4 3 6.3750 10 1 post X5 2 4.8750 11 1 pre X1 1 3.8125 12 1 pre X2 2 4.1875 13 1 pre X3 4 5.3125 14 1 pre X4 2 4.5000 15 1 pre X5 1 4.0625 16 10 fup X1 8 6.0625 17 10 fup X2 8 6.2500 18 10 fup X3 9 7.8125 19 10 fup X4 7 6.3125 20 10 fup X5 8 5.4375 21 10 post X1 6 5.0625 22 10 post X2 7 5.5625 23 10 post X3 6 6.8750 24 10 post X4 8 6.3750 25 10 post X5 8 4.8750 26 10 pre X1 4 3.8125 27 10 pre X2 4 4.1875 28 10 pre X3 5 5.3125 29 10 pre X4 3 4.5000 30 10 pre X5 4 4.0625 31 11 fup X1 5 6.0625 32 11 fup X2 6 6.2500 33 11 fup X3 8 7.8125 34 11 fup X4 6 6.3125 35 11 fup X5 5 5.4375 36 11 post X1 5 5.0625 37 11 post X2 4 5.5625 38 11 post X3 7 6.8750 39 11 post X4 5 6.3750 40 11 post X5 4 4.8750 41 11 pre X1 3 3.8125 42 11 pre X2 3 4.1875 43 11 pre X3 4 5.3125 44 11 pre X4 2 4.5000 45 11 pre X5 3 4.0625 46 12 fup X1 8 6.0625 47 12 fup X2 7 6.2500 48 12 fup X3 10 7.8125 49 12 fup X4 8 6.3125 50 12 fup X5 7 5.4375 51 12 post X1 9 5.0625 52 12 post X2 10 5.5625 53 12 post X3 11 6.8750 54 12 post X4 9 6.3750 55 12 post X5 6 4.8750 56 12 pre X1 6 3.8125 57 12 pre X2 7 4.1875 58 12 pre X3 8 5.3125 59 12 pre X4 6 4.5000 60 12 pre X5 3 4.0625 61 13 fup X1 7 6.0625 62 13 fup X2 7 6.2500 63 13 fup X3 8 7.8125 64 13 fup X4 10 6.3125 65 13 fup X5 8 5.4375 66 13 post X1 4 5.0625 67 13 post X2 6 5.5625 68 13 post X3 6 6.8750 69 13 post X4 8 6.3750 70 13 post X5 6 4.8750 71 13 pre X1 5 3.8125 72 13 pre X2 5 4.1875 73 13 pre X3 6 5.3125 74 13 pre X4 8 4.5000 75 13 pre X5 6 4.0625 76 14 fup X1 6 6.0625 77 14 fup X2 7 6.2500 78 14 fup X3 8 7.8125 79 14 fup X4 6 6.3125 80 14 fup X5 3 5.4375 81 14 post X1 5 5.0625 82 14 post X2 6 5.5625 83 14 post X3 7 6.8750 84 14 post X4 5 6.3750 85 14 post X5 2 4.8750 86 14 pre X1 2 3.8125 87 14 pre X2 2 4.1875 88 14 pre X3 3 5.3125 89 14 pre X4 1 4.5000 90 14 pre X5 2 4.0625 91 15 fup X1 7 6.0625 92 15 fup X2 7 6.2500 93 15 fup X3 8 7.8125 94 15 fup X4 6 6.3125 95 15 fup X5 7 5.4375 96 15 post X1 6 5.0625 97 15 post X2 6 5.5625 98 15 post X3 7 6.8750 99 15 post X4 9 6.3750 100 15 post X5 7 4.8750 101 15 pre X1 2 3.8125 102 15 pre X2 2 4.1875 103 15 pre X3 3 5.3125 104 15 pre X4 4 4.5000 105 15 pre X5 4 4.0625 106 16 fup X1 7 6.0625 107 16 fup X2 8 6.2500 108 16 fup X3 10 7.8125 109 16 fup X4 8 6.3125 110 16 fup X5 7 5.4375 111 16 post X1 7 5.0625 112 16 post X2 7 5.5625 113 16 post X3 8 6.8750 114 16 post X4 6 6.3750 115 16 post X5 7 4.8750 116 16 pre X1 4 3.8125 117 16 pre X2 5 4.1875 118 16 pre X3 7 5.3125 119 16 pre X4 5 4.5000 120 16 pre X5 4 4.0625 121 2 fup X1 4 6.0625 122 2 fup X2 5 6.2500 123 2 fup X3 6 7.8125 124 2 fup X4 4 6.3125 125 2 fup X5 1 5.4375 126 2 post X1 2 5.0625 127 2 post X2 2 5.5625 128 2 post X3 3 6.8750 129 2 post X4 5 6.3750 130 2 post X5 3 4.8750 131 2 pre X1 4 3.8125 132 2 pre X2 4 4.1875 133 2 pre X3 5 5.3125 134 2 pre X4 3 4.5000 135 2 pre X5 4 4.0625 136 3 fup X1 7 6.0625 137 3 fup X2 6 6.2500 138 3 fup X3 9 7.8125 139 3 fup X4 7 6.3125 140 3 fup X5 6 5.4375 141 3 post X1 4 5.0625 142 3 post X2 5 5.5625 143 3 post X3 7 6.8750 144 3 post X4 5 6.3750 145 3 post X5 4 4.8750 146 3 pre X1 5 3.8125 147 3 pre X2 6 4.1875 148 3 pre X3 5 5.3125 149 3 pre X4 7 4.5000 150 3 pre X5 7 4.0625 151 4 fup X1 4 6.0625 152 4 fup X2 4 6.2500 153 4 fup X3 5 7.8125 154 4 fup X4 3 6.3125 155 4 fup X5 4 5.4375 156 4 post X1 2 5.0625 157 4 post X2 2 5.5625 158 4 post X3 3 6.8750 159 4 post X4 5 6.3750 160 4 post X5 3 4.8750 161 4 pre X1 5 3.8125 162 4 pre X2 4 4.1875 163 4 pre X3 7 5.3125 164 4 pre X4 5 4.5000 165 4 pre X5 4 4.0625 166 5 fup X1 4 6.0625 167 5 fup X2 3 6.2500 168 5 fup X3 6 7.8125 169 5 fup X4 4 6.3125 170 5 fup X5 3 5.4375 171 5 post X1 6 5.0625 172 5 post X2 7 5.5625 173 5 post X3 8 6.8750 174 5 post X4 6 6.3750 175 5 post X5 3 4.8750 176 5 pre X1 3 3.8125 177 5 pre X2 4 4.1875 178 5 pre X3 6 5.3125 179 5 pre X4 4 4.5000 180 5 pre X5 3 4.0625 181 6 fup X1 9 6.0625 182 6 fup X2 10 6.2500 183 6 fup X3 11 7.8125 184 6 fup X4 9 6.3125 185 6 fup X5 6 5.4375 186 6 post X1 9 5.0625 187 6 post X2 9 5.5625 188 6 post X3 10 6.8750 189 6 post X4 8 6.3750 190 6 post X5 9 4.8750 191 6 pre X1 7 3.8125 192 6 pre X2 8 4.1875 193 6 pre X3 7 5.3125 194 6 pre X4 9 4.5000 195 6 pre X5 9 4.0625 196 7 fup X1 8 6.0625 197 7 fup X2 9 6.2500 198 7 fup X3 11 7.8125 199 7 fup X4 9 6.3125 200 7 fup X5 8 5.4375 201 7 post X1 7 5.0625 202 7 post X2 7 5.5625 203 7 post X3 8 6.8750 204 7 post X4 10 6.3750 205 7 post X5 8 4.8750 206 7 pre X1 5 3.8125 207 7 pre X2 5 4.1875 208 7 pre X3 6 5.3125 209 7 pre X4 4 4.5000 210 7 pre X5 5 4.0625 211 8 fup X1 6 6.0625 212 8 fup X2 6 6.2500 213 8 fup X3 7 7.8125 214 8 fup X4 5 6.3125 215 8 fup X5 6 5.4375 216 8 post X1 2 5.0625 217 8 post X2 4 5.5625 218 8 post X3 8 6.8750 219 8 post X4 6 6.3750 220 8 post X5 5 4.8750 221 8 pre X1 2 3.8125 222 8 pre X2 3 4.1875 223 8 pre X3 5 5.3125 224 8 pre X4 3 4.5000 225 8 pre X5 2 4.0625 226 9 fup X1 5 6.0625 227 9 fup X2 4 6.2500 228 9 fup X3 7 7.8125 229 9 fup X4 5 6.3125 230 9 fup X5 4 5.4375 231 9 post X1 4 5.0625 232 9 post X2 5 5.5625 233 9 post X3 6 6.8750 234 9 post X4 4 6.3750 235 9 post X5 1 4.8750 236 9 pre X1 3 3.8125 237 9 pre X2 3 4.1875 238 9 pre X3 4 5.3125 239 9 pre X4 6 4.5000 240 9 pre X5 4 4.0625 > fitted(mixed, append = TRUE) id phase hour treatment gender value .fitted 1 1 fup X1 control M 2 4.333333 2 1 fup X2 control M 3 4.666667 3 1 fup X3 control M 2 5.666667 4 1 fup X4 control M 4 5.000000 5 1 fup X5 control M 4 3.666667 6 1 post X1 control M 3 3.000000 7 1 post X2 control M 2 3.000000 8 1 post X3 control M 5 5.000000 9 1 post X4 control M 3 4.333333 10 1 post X5 control M 2 3.000000 11 1 pre X1 control M 1 3.333333 12 1 pre X2 control M 2 4.000000 13 1 pre X3 control M 4 4.666667 14 1 pre X4 control M 2 4.000000 15 1 pre X5 control M 1 4.000000 16 10 fup X1 B M 8 7.000000 17 10 fup X2 B M 8 7.000000 18 10 fup X3 B M 9 9.000000 19 10 fup X4 B M 7 7.000000 20 10 fup X5 B M 8 6.666667 21 10 post X1 B M 6 6.666667 22 10 post X2 B M 7 7.000000 23 10 post X3 B M 6 8.000000 24 10 post X4 B M 8 7.333333 25 10 post X5 B M 8 6.000000 26 10 pre X1 B M 4 4.333333 27 10 pre X2 B M 4 4.666667 28 10 pre X3 B M 5 5.666667 29 10 pre X4 B M 3 3.666667 30 10 pre X5 B M 4 3.333333 31 11 fup X1 B M 5 7.000000 32 11 fup X2 B M 6 7.000000 33 11 fup X3 B M 8 9.000000 34 11 fup X4 B M 6 7.000000 35 11 fup X5 B M 5 6.666667 36 11 post X1 B M 5 6.666667 37 11 post X2 B M 4 7.000000 38 11 post X3 B M 7 8.000000 39 11 post X4 B M 5 7.333333 40 11 post X5 B M 4 6.000000 41 11 pre X1 B M 3 4.333333 42 11 pre X2 B M 3 4.666667 43 11 pre X3 B M 4 5.666667 44 11 pre X4 B M 2 3.666667 45 11 pre X5 B M 3 3.333333 46 12 fup X1 B M 8 7.000000 47 12 fup X2 B M 7 7.000000 48 12 fup X3 B M 10 9.000000 49 12 fup X4 B M 8 7.000000 50 12 fup X5 B M 7 6.666667 51 12 post X1 B M 9 6.666667 52 12 post X2 B M 10 7.000000 53 12 post X3 B M 11 8.000000 54 12 post X4 B M 9 7.333333 55 12 post X5 B M 6 6.000000 56 12 pre X1 B M 6 4.333333 57 12 pre X2 B M 7 4.666667 58 12 pre X3 B M 8 5.666667 59 12 pre X4 B M 6 3.666667 60 12 pre X5 B M 3 3.333333 61 13 fup X1 B F 7 6.750000 62 13 fup X2 B F 7 7.250000 63 13 fup X3 B F 8 8.500000 64 13 fup X4 B F 10 7.500000 65 13 fup X5 B F 8 6.250000 66 13 post X1 B F 4 5.500000 67 13 post X2 B F 6 6.250000 68 13 post X3 B F 6 7.000000 69 13 post X4 B F 8 7.000000 70 13 post X5 B F 6 5.500000 71 13 pre X1 B F 5 3.250000 72 13 pre X2 B F 5 3.500000 73 13 pre X3 B F 6 4.750000 74 13 pre X4 B F 8 4.500000 75 13 pre X5 B F 6 4.000000 76 14 fup X1 B F 6 6.750000 77 14 fup X2 B F 7 7.250000 78 14 fup X3 B F 8 8.500000 79 14 fup X4 B F 6 7.500000 80 14 fup X5 B F 3 6.250000 81 14 post X1 B F 5 5.500000 82 14 post X2 B F 6 6.250000 83 14 post X3 B F 7 7.000000 84 14 post X4 B F 5 7.000000 85 14 post X5 B F 2 5.500000 86 14 pre X1 B F 2 3.250000 87 14 pre X2 B F 2 3.500000 88 14 pre X3 B F 3 4.750000 89 14 pre X4 B F 1 4.500000 90 14 pre X5 B F 2 4.000000 91 15 fup X1 B F 7 6.750000 92 15 fup X2 B F 7 7.250000 93 15 fup X3 B F 8 8.500000 94 15 fup X4 B F 6 7.500000 95 15 fup X5 B F 7 6.250000 96 15 post X1 B F 6 5.500000 97 15 post X2 B F 6 6.250000 98 15 post X3 B F 7 7.000000 99 15 post X4 B F 9 7.000000 100 15 post X5 B F 7 5.500000 101 15 pre X1 B F 2 3.250000 102 15 pre X2 B F 2 3.500000 103 15 pre X3 B F 3 4.750000 104 15 pre X4 B F 4 4.500000 105 15 pre X5 B F 4 4.000000 106 16 fup X1 B F 7 6.750000 107 16 fup X2 B F 8 7.250000 108 16 fup X3 B F 10 8.500000 109 16 fup X4 B F 8 7.500000 110 16 fup X5 B F 7 6.250000 111 16 post X1 B F 7 5.500000 112 16 post X2 B F 7 6.250000 113 16 post X3 B F 8 7.000000 114 16 post X4 B F 6 7.000000 115 16 post X5 B F 7 5.500000 116 16 pre X1 B F 4 3.250000 117 16 pre X2 B F 5 3.500000 118 16 pre X3 B F 7 4.750000 119 16 pre X4 B F 5 4.500000 120 16 pre X5 B F 4 4.000000 121 2 fup X1 control M 4 4.333333 122 2 fup X2 control M 5 4.666667 123 2 fup X3 control M 6 5.666667 124 2 fup X4 control M 4 5.000000 125 2 fup X5 control M 1 3.666667 126 2 post X1 control M 2 3.000000 127 2 post X2 control M 2 3.000000 128 2 post X3 control M 3 5.000000 129 2 post X4 control M 5 4.333333 130 2 post X5 control M 3 3.000000 131 2 pre X1 control M 4 3.333333 132 2 pre X2 control M 4 4.000000 133 2 pre X3 control M 5 4.666667 134 2 pre X4 control M 3 4.000000 135 2 pre X5 control M 4 4.000000 136 3 fup X1 control M 7 4.333333 137 3 fup X2 control M 6 4.666667 138 3 fup X3 control M 9 5.666667 139 3 fup X4 control M 7 5.000000 140 3 fup X5 control M 6 3.666667 141 3 post X1 control M 4 3.000000 142 3 post X2 control M 5 3.000000 143 3 post X3 control M 7 5.000000 144 3 post X4 control M 5 4.333333 145 3 post X5 control M 4 3.000000 146 3 pre X1 control M 5 3.333333 147 3 pre X2 control M 6 4.000000 148 3 pre X3 control M 5 4.666667 149 3 pre X4 control M 7 4.000000 150 3 pre X5 control M 7 4.000000 151 4 fup X1 control F 4 4.000000 152 4 fup X2 control F 4 3.500000 153 4 fup X3 control F 5 5.500000 154 4 fup X4 control F 3 3.500000 155 4 fup X5 control F 4 3.500000 156 4 post X1 control F 2 4.000000 157 4 post X2 control F 2 4.500000 158 4 post X3 control F 3 5.500000 159 4 post X4 control F 5 5.500000 160 4 post X5 control F 3 3.000000 161 4 pre X1 control F 5 4.000000 162 4 pre X2 control F 4 4.000000 163 4 pre X3 control F 7 6.500000 164 4 pre X4 control F 5 4.500000 165 4 pre X5 control F 4 3.500000 166 5 fup X1 control F 4 4.000000 167 5 fup X2 control F 3 3.500000 168 5 fup X3 control F 6 5.500000 169 5 fup X4 control F 4 3.500000 170 5 fup X5 control F 3 3.500000 171 5 post X1 control F 6 4.000000 172 5 post X2 control F 7 4.500000 173 5 post X3 control F 8 5.500000 174 5 post X4 control F 6 5.500000 175 5 post X5 control F 3 3.000000 176 5 pre X1 control F 3 4.000000 177 5 pre X2 control F 4 4.000000 178 5 pre X3 control F 6 6.500000 179 5 pre X4 control F 4 4.500000 180 5 pre X5 control F 3 3.500000 181 6 fup X1 A M 9 8.500000 182 6 fup X2 A M 10 9.500000 183 6 fup X3 A M 11 11.000000 184 6 fup X4 A M 9 9.000000 185 6 fup X5 A M 6 7.000000 186 6 post X1 A M 9 8.000000 187 6 post X2 A M 9 8.000000 188 6 post X3 A M 10 9.000000 189 6 post X4 A M 8 9.000000 190 6 post X5 A M 9 8.500000 191 6 pre X1 A M 7 6.000000 192 6 pre X2 A M 8 6.500000 193 6 pre X3 A M 7 6.500000 194 6 pre X4 A M 9 6.500000 195 6 pre X5 A M 9 7.000000 196 7 fup X1 A M 8 8.500000 197 7 fup X2 A M 9 9.500000 198 7 fup X3 A M 11 11.000000 199 7 fup X4 A M 9 9.000000 200 7 fup X5 A M 8 7.000000 201 7 post X1 A M 7 8.000000 202 7 post X2 A M 7 8.000000 203 7 post X3 A M 8 9.000000 204 7 post X4 A M 10 9.000000 205 7 post X5 A M 8 8.500000 206 7 pre X1 A M 5 6.000000 207 7 pre X2 A M 5 6.500000 208 7 pre X3 A M 6 6.500000 209 7 pre X4 A M 4 6.500000 210 7 pre X5 A M 5 7.000000 211 8 fup X1 A F 6 5.500000 212 8 fup X2 A F 6 5.000000 213 8 fup X3 A F 7 7.000000 214 8 fup X4 A F 5 5.000000 215 8 fup X5 A F 6 5.000000 216 8 post X1 A F 2 3.000000 217 8 post X2 A F 4 4.500000 218 8 post X3 A F 8 7.000000 219 8 post X4 A F 6 5.000000 220 8 post X5 A F 5 3.000000 221 8 pre X1 A F 2 2.500000 222 8 pre X2 A F 3 3.000000 223 8 pre X3 A F 5 4.500000 224 8 pre X4 A F 3 4.500000 225 8 pre X5 A F 2 3.000000 226 9 fup X1 A F 5 5.500000 227 9 fup X2 A F 4 5.000000 228 9 fup X3 A F 7 7.000000 229 9 fup X4 A F 5 5.000000 230 9 fup X5 A F 4 5.000000 231 9 post X1 A F 4 3.000000 232 9 post X2 A F 5 4.500000 233 9 post X3 A F 6 7.000000 234 9 post X4 A F 4 5.000000 235 9 post X5 A F 1 3.000000 236 9 pre X1 A F 3 2.500000 237 9 pre X2 A F 3 3.000000 238 9 pre X3 A F 4 4.500000 239 9 pre X4 A F 6 4.500000 240 9 pre X5 A F 4 3.000000 > fitted(between, append = TRUE) id treatment gender value .fitted 1 1 control M 2.666667 4.111111 2 2 control M 3.666667 4.111111 3 3 control M 6.000000 4.111111 4 4 control F 4.000000 4.333333 5 5 control F 4.666667 4.333333 6 6 A M 8.666667 8.000000 7 7 A M 7.333333 8.000000 8 8 A F 4.666667 4.500000 9 9 A F 4.333333 4.500000 10 10 B M 6.333333 6.222222 11 11 B M 4.666667 6.222222 12 12 B M 7.666667 6.222222 13 13 B F 6.666667 5.833333 14 14 B F 4.333333 5.833333 15 15 B F 5.666667 5.833333 16 16 B F 6.666667 5.833333 > > ## No message: > fitted(between2) 1 2 3 4 5 6 7 8 4.111111 4.111111 4.111111 4.333333 4.333333 8.000000 8.000000 4.500000 9 10 11 12 13 14 15 16 4.500000 6.222222 6.222222 6.222222 5.833333 5.833333 5.833333 5.833333 > fitted(within2) 1 2 3 4 5 6 7 8 9 10 11 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 12 13 14 15 16 17 18 19 20 21 22 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 23 24 25 26 27 28 29 30 31 32 33 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 34 35 36 37 38 39 40 41 42 43 44 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 45 46 47 48 49 50 51 52 53 54 55 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 56 57 58 59 60 61 62 63 64 65 66 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 67 68 69 70 71 72 73 74 75 76 77 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 78 79 80 81 82 83 84 85 86 87 88 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 89 90 91 92 93 94 95 96 97 98 99 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 100 101 102 103 104 105 106 107 108 109 110 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 111 112 113 114 115 116 117 118 119 120 121 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 122 123 124 125 126 127 128 129 130 131 132 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 133 134 135 136 137 138 139 140 141 142 143 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 144 145 146 147 148 149 150 151 152 153 154 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 155 156 157 158 159 160 161 162 163 164 165 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 166 167 168 169 170 171 172 173 174 175 176 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 177 178 179 180 181 182 183 184 185 186 187 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 188 189 190 191 192 193 194 195 196 197 198 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 199 200 201 202 203 204 205 206 207 208 209 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 210 211 212 213 214 215 216 217 218 219 220 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 5.5625 6.8750 6.3750 4.8750 221 222 223 224 225 226 227 228 229 230 231 3.8125 4.1875 5.3125 4.5000 4.0625 6.0625 6.2500 7.8125 6.3125 5.4375 5.0625 232 233 234 235 236 237 238 239 240 5.5625 6.8750 6.3750 4.8750 3.8125 4.1875 5.3125 4.5000 4.0625 > > #### residuals() and fitted() methods can be used for plotting > ### requires package ggResidpanel > if (require("ggResidpanel")) { + resid_auxpanel(residuals = residuals(mixed), predicted = fitted(mixed)) + + ## Not run: + ##D ## suppress Messages: + ##D suppressMessages( + ##D resid_auxpanel(residuals = residuals(mixed), predicted = fitted(mixed)) + ##D ) + ##D + ## End(Not run) + } Loading required package: ggResidpanel Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data. residuals(..., append = TRUE) will return data and residuals. Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data. fitted(..., append = TRUE) will return data and fitted values. > > > > cleanEx() detaching ‘package:ggResidpanel’ > nameEx("round_ps") > ### * round_ps > > flush(stderr()); flush(stdout()) > > ### Name: round_ps > ### Title: Helper functions for rounding p-values > ### Aliases: round_ps round_ps_apa > > ### ** Examples > > x <- runif(10) > y <- runif(10, 0, .01) > > round_ps(x) [1] ".27" ".37" ".57" ".91" ".20" ".90" ".94" ".66" ".63" ".06" > round_ps_apa(x) [1] ".266" ".372" ".573" ".908" ".202" ".898" ".945" ".661" ".629" ".062" > > round_ps(y) [1] ".002" ".002" ".007" ".004" ".008" ".005" ".007" ".010" ".004" ".008" > round_ps_apa(y) [1] ".002" ".002" ".007" ".004" ".008" ".005" ".007" ".010" ".004" ".008" > > round_ps(0.0000000099) [1] "<.0001" > round_ps_apa(0.0000000099) [1] "<.001" > > > > > cleanEx() > nameEx("sk2011.1") > ### * sk2011.1 > > flush(stderr()); flush(stdout()) > > ### Name: sk2011.1 > ### Title: Data from Singmann & Klauer (2011, Experiment 1) > ### Aliases: sk2011.1 > ### Keywords: dataset > > ### ** Examples > > data(sk2011.1) > > # Table 1 (p. 264): > aov_ez("id", "response", sk2011.1[ sk2011.1$what == "affirmation",], + within = c("inference", "type"), between = "instruction", + anova_table=(es = "pes")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: instruction Anova Table (Type 3 tests) Response: response Effect df MSE F pes p.value 1 instruction 1, 38 1072.42 0.13 .003 .723 2 inference 1, 38 1007.21 13.01 *** .255 <.001 3 instruction:inference 1, 38 1007.21 12.44 ** .247 .001 4 type 1, 38 187.90 0.06 .002 .805 5 instruction:type 1, 38 187.90 3.09 + .075 .087 6 inference:type 1, 38 498.48 29.62 *** .438 <.001 7 instruction:inference:type 1, 38 498.48 10.73 ** .220 .002 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > aov_ez("id", "response", sk2011.1[ sk2011.1$what == "denial",], + within = c("inference", "type"), between = "instruction", + anova_table=(es = "pes")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: instruction Anova Table (Type 3 tests) Response: response Effect df MSE F pes p.value 1 instruction 1, 38 1615.52 1.37 .035 .250 2 inference 1, 38 884.57 1.78 .045 .190 3 instruction:inference 1, 38 884.57 1.19 .030 .281 4 type 1, 38 338.60 0.13 .003 .720 5 instruction:type 1, 38 338.60 3.72 + .089 .061 6 inference:type 1, 38 174.88 18.99 *** .333 <.001 7 instruction:inference:type 1, 38 174.88 4.13 * .098 .049 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 > > > > > > cleanEx() > nameEx("sk2011.2") > ### * sk2011.2 > > flush(stderr()); flush(stdout()) > > ### Name: sk2011.2 > ### Title: Data from Singmann & Klauer (2011, Experiment 2) > ### Aliases: sk2011.2 > ### Keywords: dataset > > ### ** Examples > > data("sk2011.2") > > ## remove excluded participants: > > sk2_final <- droplevels(sk2011.2[!(sk2011.2$id %in% c(7, 8, 9, 12, 16, 17, 24, 30)),]) > str(sk2_final) 'data.frame': 1980 obs. of 8 variables: $ id : Factor w/ 55 levels "1","2","3","4",..: 1 3 6 8 10 13 15 17 20 21 ... $ instruction: Factor w/ 2 levels "deductive","inductive": 1 1 1 1 1 1 1 1 1 1 ... $ inference : Factor w/ 4 levels "MP","MT","AC",..: 1 1 1 1 1 1 1 1 1 1 ... $ validity : Factor w/ 2 levels "valid","invalid": 1 1 1 1 1 1 1 1 1 1 ... $ what : Factor w/ 2 levels "affirmation",..: 1 1 1 1 1 1 1 1 1 1 ... $ type : Factor w/ 3 levels "prological","neutral",..: 1 1 1 1 1 1 1 1 1 1 ... $ response : int 100 100 100 100 99 99 99 100 100 100 ... $ content : Factor w/ 9 levels "C1","C2","C3",..: 1 1 1 1 1 1 1 1 1 1 ... > > ## Table 2 (inference = problem): > aov_ez("id", "response", sk2_final[sk2_final$what == "affirmation",], + between = "instruction", within = c("inference", "type"), + anova_table=list(es = "pes")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: instruction Anova Table (Type 3 tests) Response: response Effect df MSE F pes p.value 1 instruction 1, 53 894.54 0.01 <.001 .943 2 inference 1, 53 978.22 87.13 *** .622 <.001 3 instruction:inference 1, 53 978.22 19.78 *** .272 <.001 4 type 1.87, 99.05 180.76 29.32 *** .356 <.001 5 instruction:type 1.87, 99.05 180.76 3.83 * .067 .028 6 inference:type 1.72, 91.18 267.53 63.93 *** .547 <.001 7 instruction:inference:type 1.72, 91.18 267.53 4.84 * .084 .014 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > aov_ez("id", "response", sk2_final[sk2_final$what == "denial",], + between = "instruction", within = c("inference", "type"), + anova_table=list(es = "pes")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: instruction Anova Table (Type 3 tests) Response: response Effect df MSE F pes p.value 1 instruction 1, 53 852.26 1.33 .024 .255 2 inference 1, 53 1300.62 4.67 * .081 .035 3 instruction:inference 1, 53 1300.62 0.47 .009 .498 4 type 1.98, 104.74 297.76 22.01 *** .293 <.001 5 instruction:type 1.98, 104.74 297.76 1.51 .028 .226 6 inference:type 1.95, 103.16 248.62 59.78 *** .530 <.001 7 instruction:inference:type 1.95, 103.16 248.62 3.09 + .055 .051 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 Sphericity correction method: GG > > # Recreate Figure 4 (corrected version): > > sk2_aff <- droplevels(sk2_final[sk2_final$what == "affirmation",]) > sk2_aff$type2 <- factor(sk2_aff$inference:sk2_aff$type, levels = c("MP:prological", + "MP:neutral", "MP:counterlogical", "AC:counterlogical", + "AC:neutral", "AC:prological")) > a1_b <- aov_ez("id", "response", sk2_aff, + between = "instruction", within = c("type2")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: instruction > > sk2_den <- droplevels(sk2_final[sk2_final$what == "denial",]) > sk2_den$type2 <- factor(sk2_den$inference:sk2_den$type, levels = c("MT:prological", + "MT:neutral", "MT:counterlogical", "DA:counterlogical", + "DA:neutral","DA:prological")) > a2_b <- aov_ez("id", "response", sk2_den, + between = "instruction", within = c("type2")) Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`. To turn off this warning, pass `fun_aggregate = mean` explicitly. Contrasts set to contr.sum for the following variables: instruction > > if (requireNamespace("emmeans") && requireNamespace("ggplot2")) { + afex_plot(a1_b,"type2", "instruction") + + ggplot2::coord_cartesian(ylim = c(0, 100)) + afex_plot(a2_b,"type2", "instruction") + + ggplot2::coord_cartesian(ylim = c(0, 100)) + } Warning: Panel(s) show a mixed within-between-design. Error bars do not allow comparisons across all means. Suppress error bars with: error = "none" Warning: Panel(s) show a mixed within-between-design. Error bars do not allow comparisons across all means. Suppress error bars with: error = "none" > > > > ### *