## ----echo=FALSE--------------------------------------------------------------- req_suggested_packages <- c("emmeans", "multcomp", "dplyr", "tidyr","ggplot2") pcheck <- lapply(req_suggested_packages, requireNamespace, quietly = TRUE) if (any(!unlist(pcheck))) { message("Required package(s) for this vignette are not available/installed and code will not be executed.") knitr::opts_chunk$set(eval = FALSE) } ## ----set-options, echo=FALSE, cache=FALSE----------------------------------------------- op <- options(width = 90, dplyr.summarise.inform = FALSE) knitr::opts_chunk$set(dpi=72) ## ----echo=FALSE------------------------------------------------------------------------- load(system.file("extdata/", "output_mixed_vignette.rda", package = "afex")) ## ----message=FALSE, warning=FALSE------------------------------------------------------- library("afex") # needed for mixed() and attaches lme4 automatically. library("emmeans") # emmeans is needed for follow-up tests library("multcomp") # for advanced control for multiple testing/Type 1 errors. library("dplyr") # for working with data frames library("tidyr") # for transforming data frames from wide to long and the other way round. library("ggplot2") # for plots theme_set(theme_bw(base_size = 15) + theme(legend.position="bottom", panel.grid.major.x = element_blank())) data("fhch2010") # load fhch <- droplevels(fhch2010[ fhch2010$correct,]) # remove errors str(fhch2010) # structure of the data ## --------------------------------------------------------------------------------------- ## are all participants in only one task? fhch %>% group_by(id) %>% summarise(task = n_distinct(task)) %>% as.data.frame() %>% {.$task == 1} %>% all() ## participants per condition: fhch %>% group_by(id) %>% summarise(task = first(task)) %>% ungroup() %>% group_by(task) %>% summarise(n = n()) ## number of different items per participant: fhch %>% group_by(id, stimulus) %>% summarise(items = n_distinct(item)) %>% ungroup() %>% group_by(stimulus) %>% summarise(min = min(items), max = max(items), mean = mean(items)) ## ----fig.width=7, fig.height=4---------------------------------------------------------- fhch_long <- fhch %>% pivot_longer(cols = c(rt, log_rt), names_to = "rt_type", values_to = "rt") ggplot(fhch_long, aes(rt)) + geom_histogram(bins = 100) + facet_wrap(vars(rt_type), scales = "free_x") ## ----fig.width=7, fig.height=6, message=FALSE, out.width="90%"-------------------------- agg_p <- fhch %>% group_by(id, task, stimulus, density, frequency) %>% summarise(mean = mean(log_rt)) %>% ungroup() ggplot(agg_p, aes(x = interaction(density,frequency), y = mean)) + ggbeeswarm::geom_quasirandom(alpha = 0.5) + geom_boxplot(fill = "transparent") + stat_summary(colour = "red") + facet_grid(cols = vars(task), rows = vars(stimulus)) ## ----fig.width=7, fig.height=6, message=FALSE, out.width="90%"-------------------------- agg_i <- fhch %>% group_by(item, task, stimulus, density, frequency) %>% summarise(mean = mean(log_rt)) %>% ungroup() ggplot(agg_i, aes(x = interaction(density,frequency), y = mean)) + ggbeeswarm::geom_quasirandom(alpha = 0.3) + geom_boxplot(fill = "transparent") + stat_summary(colour = "red") + facet_grid(cols = vars(task), rows = vars(stimulus)) ## ----eval = FALSE----------------------------------------------------------------------- # m1s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus*density*frequency|id)+ # (task|item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6))) ## ----echo=FALSE------------------------------------------------------------------------- message("boundary (singular) fit: see ?isSingular") ## ----eval = FALSE----------------------------------------------------------------------- # m2s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus*density*frequency|id)+ # (task||item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----echo=FALSE------------------------------------------------------------------------- message("boundary (singular) fit: see ?isSingular") ## ----eval = FALSE----------------------------------------------------------------------- # m3s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus*density*frequency||id)+ # (task|item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----echo=FALSE------------------------------------------------------------------------- message("boundary (singular) fit: see ?isSingular") ## ----eval = FALSE----------------------------------------------------------------------- # m4s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus*density*frequency||id)+ # (task||item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----echo=FALSE------------------------------------------------------------------------- message("boundary (singular) fit: see ?isSingular") ## ----eval=FALSE------------------------------------------------------------------------- # summary(m4s)$varcor ## ----echo=FALSE------------------------------------------------------------------------- cat(outp_m4s_vc$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # m5s <- mixed(log_rt ~ task*stimulus*density*frequency + # ((stimulus+density+frequency)^2||id)+ # (task||item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----echo=FALSE------------------------------------------------------------------------- message("boundary (singular) fit: see ?isSingular") ## ----eval=FALSE------------------------------------------------------------------------- # summary(m5s)$varcor ## ----echo=FALSE------------------------------------------------------------------------- cat(outp_m5s_vc$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # m6s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus+density+frequency||id)+ # (task||item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), # expand_re = TRUE) ## ----echo=FALSE------------------------------------------------------------------------- message("boundary (singular) fit: see ?isSingular") ## ----eval=FALSE------------------------------------------------------------------------- # summary(m6s)$varcor ## ----echo=FALSE------------------------------------------------------------------------- cat(outp_m6s_vc$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # m7s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus+frequency||id)+ # (task||item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----eval=FALSE------------------------------------------------------------------------- # m8s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus+frequency|id)+ # (task||item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----echo=FALSE------------------------------------------------------------------------- warning(fit_m8s$warnings, call. = FALSE) ## ----eval=FALSE------------------------------------------------------------------------- # m9s <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus+frequency||id)+ # (task|item), fhch, # control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) ## ----eval=FALSE------------------------------------------------------------------------- # left_join(nice(m1s), nice(m9s), by = "Effect", # suffix = c("_full", "_final")) ## ----echo=FALSE------------------------------------------------------------------------- cat(outp_comb_anova$output, sep = "\n") ## ----eval = FALSE----------------------------------------------------------------------- # m9lrt <- mixed(log_rt ~ task*stimulus*density*frequency + # (stimulus+frequency||id)+ # (task|item), fhch, method = "LRT", # control = lmerControl(optCtrl = list(maxfun = 1e6)), # expand_re = TRUE) # m9lrt ## ----echo=FALSE------------------------------------------------------------------------- cat(outp_m9lrt$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # emm_options(lmer.df = "asymptotic") # also possible: 'satterthwaite', 'kenward-roger' # emm_i1 <- emmeans(m9s, "frequency", by = c("stimulus", "task")) # emm_i1 ## ----echo=FALSE------------------------------------------------------------------------- message("NOTE: Results may be misleading due to involvement in interactions") cat(emm_o1$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # update(pairs(emm_i1), by = NULL, adjust = "holm") ## ----echo=FALSE------------------------------------------------------------------------- message("NOTE: Results may be misleading due to involvement in interactions") cat(emm_o2$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # summary(as.glht(update(pairs(emm_i1), by = NULL)), test = adjusted("free")) ## ----echo=FALSE------------------------------------------------------------------------- cat(emm_o3$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # emm_i1b <- update(emm_i1, tran = "log", type = "response", by = NULL) # emm_i1b ## ----echo=FALSE------------------------------------------------------------------------- cat(emm_o4$output, sep = "\n") ## ----echo=FALSE------------------------------------------------------------------------- load(system.file("extdata/", "output_afex_plot_mixed_vignette_model.rda", package = "afex")) emm_options(lmer.df = "asymptotic") ## ----eval=TRUE, fig.width=5, fig.height=4, out.width="90%"------------------------------ afex_plot(m9s, "frequency", "stimulus", "task", id = "id", data_geom = ggbeeswarm::geom_quasirandom, data_arg = list( dodge.width = 0.5, ## needs to be same as dodge cex = 0.8, width = 0.1, color = "darkgrey")) ## ----eval=FALSE------------------------------------------------------------------------- # joint_tests(m9s, by = c("stimulus", "task")) ## ----echo=FALSE------------------------------------------------------------------------- cat(emm_o5$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # emm_i2 <- emmeans(m2s, c("density", "frequency"), by = c("stimulus", "task")) # emm_i2 ## ----echo=FALSE------------------------------------------------------------------------- cat(emm_o6$output, sep = "\n") ## ----eval=FALSE------------------------------------------------------------------------- # # desired contrats: # des_c <- list( # ll_hl = c(1, -1, 0, 0), # lh_hh = c(0, 0, 1, -1) # ) # update(contrast(emm_i2, des_c), by = NULL, adjust = "holm") ## ----echo=FALSE------------------------------------------------------------------------- cat(emm_o7$output, sep = "\n") ## ----include=FALSE------------------------------------------------------------ options(op)