tests/testthat/test_apa_print_emm_lsm.R
In crsh/papaja: Prepare American Psychological Association Journal Articles with R Markdown
context("apa_print() for emmeans/lsmeans")
test_that(
"ANOVA"
, {
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
))
# Main effect ------------------------------------------------------
tw_me_emm <- emmeans::emmeans(tw_rm, ~ Valence)
tw_me_emm_output <- expect_warning(
apa_print(tw_me_emm, conf.l = .95)
, regexp = "Using argument 'conf.level' in calls to 'apa_print()' is deprecated. Please use 'conf.int' instead."
, fixed = TRUE
)
expect_apa_results(
tw_me_emm_output
, labels = list(
Valence = "Valence"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = levels(tw_rm_data$Valence)
)
expect_apa_term(
tw_me_emm_output
, term = "Neg"
, estimate = "$M = 11.00$, 95\\% CI $[6.34, 15.66]$"
, statistic = "$t(4.63) = 6.21$, $p = .002$"
)
# Alternative calls
tw_me_lsm <- emmeans::lsmeans(tw_rm, ~ Valence)
tw_me_lsm_output <- apa_print(tw_me_lsm)
tw_me_lsm_output2 <- apa_print(
summary(tw_me_lsm, infer = TRUE)
, est_name = "M"
)
tw_me_emm_output2 <- apa_print(
summary(tw_me_emm, infer = TRUE)
, est_name = "M"
)
expect_identical(tw_me_emm_output, tw_me_emm_output2)
expect_identical(tw_me_lsm_output, tw_me_lsm_output2)
expect_identical(tw_me_lsm_output, tw_me_emm_output2)
expect_warning(
tw_me_emm_p_output <- apa_print(
summary(tw_me_emm, infer = c(FALSE, TRUE))
, est_name = "M"
)
)
expect_identical(
tw_me_lsm_output$table[, c("Valence", "estimate", "statistic", "df", "p.value")]
, tw_me_emm_p_output$table
)
tw_me_emm_ci_output <- apa_print(
summary(tw_me_emm, infer = c(TRUE, FALSE))
, est_name = "M"
)
expect_apa_term(
tw_me_emm_ci_output
, term = "Neg"
, estimate = "$M = 11.00$, 95\\% CI $[6.34, 15.66]$"
, statistic = NULL
)
expect_identical(
tw_me_lsm_output$table[, c("Valence", "estimate", "conf.int")]
, tw_me_emm_ci_output$table
)
expect_apa_results(
apa_print(summary(tw_me_lsm, infer = c(FALSE, FALSE)))
, labels = list(
Valence = "Valence"
, estimate = "$\\hat{\\theta}$"
)
, term_names = levels(tw_rm_data$Valence)
)
tw_me_emm_output3 <- apa_print(
summary(tw_me_emm, infer = TRUE, side = ">")
)
# https://github.com/crsh/papaja/issues/456#issuecomment-901653372
expect_apa_term(
tw_me_emm_output3
, term = "Neg"
, estimate = "$\\hat{\\theta} = 11.00$, 95\\% CI $[7.37, \\infty]$"
, statistic = "$t(4.63) = 6.21$, $p = .001$"
)
# Interaction ------------------------------------------------------
tw_int_emm <- emmeans::emmeans(tw_rm, ~ Valence * Task)
tw_int_emm_output <- apa_print(tw_int_emm)
term_names <- apply(expand.grid(levels(tw_rm_data$Valence), levels(tw_rm_data$Task)), 1, paste, collapse = "_")
expect_apa_results(
tw_int_emm_output
, labels = list(
Task = "Task"
, Valence = "Valence"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = term_names
)
expect_apa_term(
tw_int_emm_output
, term = "Neg_Cued"
, estimate = "$M = 11.80$, 95\\% CI $[7.17, 16.43]$"
, statistic = "$t(5.52) = 6.37$, $p < .001$"
)
# Alternative calls
tw_int_lsm <- emmeans::lsmeans(tw_rm, ~ Valence * Task)
tw_int_lsm_output <- apa_print(tw_int_lsm)
tw_int_lsm_output2 <- apa_print(
summary(tw_int_lsm, infer = TRUE)
, est_name = "M"
)
tw_int_emm_output2 <- apa_print(
summary(tw_int_emm, infer = TRUE)
, est_name = "M"
)
expect_identical(tw_int_emm_output, tw_int_emm_output2)
expect_identical(tw_int_lsm_output, tw_int_lsm_output2)
expect_identical(tw_int_lsm_output, tw_int_emm_output2)
# Simple effects ---------------------------------------------------
tw_se_emm <- emmeans::emmeans(tw_rm, ~ Valence | Task)
tw_se_emm_output <- apa_print(tw_se_emm)
expect_apa_results(
tw_se_emm_output
, labels = list(
Task = "Task"
, Valence = "Valence"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = term_names
)
# Sort
tw_int_emm_output$estimate <- tw_int_emm_output$estimate[names(tw_se_emm_output$estimate)]
tw_int_emm_output$statistic <- tw_int_emm_output$statistic[names(tw_int_emm_output$estimate)]
tw_int_emm_output$full_result <- tw_int_emm_output$full_result[names(tw_int_emm_output$estimate)]
tw_int_emm_output$table <- tw_int_emm_output$table[rownames(tw_se_emm_output$table), ]
expect_identical(tw_se_emm_output, tw_int_emm_output)
# Complex output (four factors) ------------------------------------
data(obk.long, package = "afex")
fw_mixed <- suppressWarnings(afex::aov_ez(
id = "id"
, dv = "value"
, data = obk.long
, between = c("treatment", "gender")
, within = c("phase", "hour")
, observed = "gender"
))
fw_mixed_emm <- emmeans::emmeans(fw_mixed, ~ phase * hour | treatment * gender)
fw_mixed_emm_output <- apa_print(fw_mixed_emm)
expect_apa_results(
fw_mixed_emm_output
, labels = list(
gender = "gender"
, treatment = "treatment"
, hour = "hour"
, phase = "phase"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
)
# table ------------------------------------------------------------
expect_identical(
object = fw_mixed_emm_output$table$gender
, expected = structure(
c(
"F"
, rep("", nrow(fw_mixed_emm_output$table) / 2 - 1)
, "M"
, rep("", nrow(fw_mixed_emm_output$table) / 2 - 1)
)
, class = c("tiny_labelled", "character")
, label = "gender"
)
)
expect_identical(
object = fw_mixed_emm_output$table$treatment
, expected = structure(
rep(
c(
"Control"
, rep("", nrow(fw_mixed_emm_output$table) / 6 - 1)
, "A"
, rep("", nrow(fw_mixed_emm_output$table) / 6 - 1)
, "B"
, rep("", nrow(fw_mixed_emm_output$table) / 6 - 1)
)
, 2
)
, class = c("tiny_labelled", "character")
, label = "treatment"
)
)
expect_identical(
object = fw_mixed_emm_output$table$hour
, expected = structure(
rep(
c(
"X1"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X2"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X3"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X4"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X5"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
)
, 6
)
, class = c("tiny_labelled", "character")
, label = "hour"
)
)
expect_identical(
object = fw_mixed_emm_output$table$phase
, expected = structure(
rep(
c("Fup", "Post", "Pre")
, 30
)
, class = c("tiny_labelled", "character")
, label = "phase"
)
)
# Alternative calls
fw_mixed_lsm <- emmeans::lsmeans(fw_mixed, ~ phase * hour | treatment * gender)
fw_mixed_lsm_output <- apa_print(fw_mixed_lsm)
fw_mixed_lsm_output2 <- apa_print(
summary(fw_mixed_lsm, infer = TRUE)
, est_name = "M"
)
fw_mixed_emm_output2 <- apa_print(
summary(fw_mixed_emm, infer = TRUE)
, est_name = "M"
)
expect_identical(fw_mixed_lsm_output, fw_mixed_lsm_output2)
expect_identical(fw_mixed_emm_output, fw_mixed_emm_output2)
expect_identical(fw_mixed_lsm_output, fw_mixed_emm_output2)
# Pairs ------------------------------------------------------------
tw_me_pairs_emm <- pairs(tw_me_emm)
tw_me_pairs_emm_output <- apa_print(tw_me_pairs_emm)
expect_apa_results(
tw_me_pairs_emm_output
, labels = list(
contrast = "Contrast"
, estimate = "$\\Delta M$"
, conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Tukey(3)}$"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Tukey(3)}$"
)
, term_names = c("Neg_Neu", "Neg_Pos", "Neu_Pos")
)
expect_apa_term(
tw_me_pairs_emm_output
, term = "Neg_Neu"
, estimate = "$\\Delta M = -1.10$, $95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Tukey(3)}$ $[-3.44, 1.24]$"
, statistic = "$t(8) = -1.34$, $p_\\mathrm{\\scriptstyle Tukey(3)} = .413$"
)
## Custom contrast names
tw_me_pairs_emm_output <- apa_print(
tw_me_pairs_emm
, contrast_names = letters[1:3]
)
expect_identical(
object = tw_me_pairs_emm_output$table$contrast
, expected = structure(
letters[1:3]
, label = "Contrast"
, class = c("tiny_labelled", "character")
)
)
# Custom contrasts
tw_se_contrast_emm <- emmeans::contrast(
tw_se_emm
, method = list("Positive - Negative" = c(1, 0, -1))
)
tw_se_contrast_emm_output <- apa_print(tw_se_contrast_emm)
tw_se_contrast_emm_output2 <- apa_print(
summary(tw_se_contrast_emm, infer = TRUE)
, est_name = "\\Delta M"
)
expect_identical(
tw_se_contrast_emm_output
, tw_se_contrast_emm_output2
)
# Simple contrasts
simple_pairs <- pairs(
emmeans::emmeans(fw_mixed, ~ gender | treatment)
)
tw_between_emm <- emmeans::emmeans(fw_mixed, ~ treatment * gender)
simple_pairs2 <- pairs(
tw_between_emm
, simple = "gender"
)
simple_pairs_output <- apa_print(simple_pairs)
simple_pairs2_output <- apa_print(simple_pairs2)
expect_apa_results(
simple_pairs_output
, term_names = c("Control_F_M", "A_F_M", "B_F_M")
)
expect_identical(simple_pairs_output, simple_pairs2_output)
simple_contrasts <- emmeans::contrast(
tw_between_emm
, "consec"
, simple = "each"
, combine = TRUE
, adjust = "none"
)
simple_contrasts_output <- apa_print(simple_contrasts)
simple_contrasts2 <- emmeans::contrast(
tw_between_emm
, "consec"
, simple = "each"
, adjust = "none"
)
simple_contrasts2_output1 <- apa_print(simple_contrasts2$`simple contrasts for gender`)
simple_contrasts2_output2 <- apa_print(simple_contrasts2$`simple contrasts for treatment`)
expect_identical(
simple_contrasts_output$full_result[grepl("M_F$", names(simple_contrasts_output$full_result))]
, simple_contrasts2_output1$full_result
)
expect_identical(
simple_contrasts_output$full_result[grepl("^(F|M)_", names(simple_contrasts_output$full_result))]
, simple_contrasts2_output2$full_result
)
}
)
test_that(
"Joint tests"
, {
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
# Joint tests
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
, anova_table = list(correction = "none")
))
tw_rm_output <- suppressWarnings(apa_print(tw_rm))
tw_rm_emm <- emmeans::emmeans(tw_rm$aov, ~ Task * Valence)
## All terms
emm_aov <- emmeans::joint_tests(tw_rm_emm)
emm_aov_output <- apa_print(emm_aov)
expect_apa_results(
emm_aov_output
, labels = list(
term = "Effect"
, statistic = "$F$"
, df = "$\\mathit{df}$"
, df.residual = "$\\mathit{df}_{\\mathrm{res}}$"
, p.value = "$p$"
)
, term_names = names(tw_rm_output$estimate)
)
expect_apa_term(
emm_aov_output
, term = "Task"
, estimate = NULL
, statistic = tw_rm_output$statistic$Task
)
## Split by
emm_split_aov <- emmeans::joint_tests(tw_rm_emm, by = "Task")
emm_split_aov_output <- apa_print(emm_split_aov)
expect_apa_results(
emm_split_aov_output
, labels = list(
Task = "Task"
, term = "Effect"
, statistic = "$F$"
, df = "$\\mathit{df}$"
, df.residual = "$\\mathit{df}_{\\mathrm{res}}$"
, p.value = "$p$"
)
, term_names = paste("Valence", levels(tw_rm_data$Task), sep = "_")
, table_terms = rep("Valence", 2) # Levels split by are in "Task" column
)
expect_apa_term(
emm_split_aov_output
, term = "Valence_Cued"
, estimate = NULL
, statistic = "$F(2, 15.58) = 1.46$, $p = .263$"
)
# Ensure proper sorting of terms
load("data/mixed_data.rdata")
unsorted_aov <- afex::aov_4(formula = Recall ~ Gender * Dosage * (Task * Valence |Subject), data = mixed_data, fun_aggregate = mean)
unsorted_emm <- emmeans::joint_tests(unsorted_aov, by = "Gender")
apa_out <- apa_print(unsorted_emm)
expect_apa_results(
apa_out
, labels = list(
Gender = "Gender"
, term = "Effect"
, statistic = "$F$"
, df = "$\\mathit{df}$"
, df.residual = "$\\mathit{df}_{\\mathrm{res}}$"
, p.value = "$p$"
)
, term_names = papaja:::sanitize_terms(paste(unlabel(gsub(apa_out$table$term, pattern = " $\\times$ ", replacement = "_", fixed = TRUE)), apa_out$table$Gender, sep = "_"))
, table_terms = beautify_terms(data.frame(unsorted_emm)$model.term)
)
}
)
test_that(
"Estimate name guessing"
, {
# ANOVA
library("emmeans")
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
))
ow_me_lsm <- lsmeans(tw_rm, ~ Valence)
expect_identical(est_name_from_call(ow_me_lsm), "M")
tw_int_lsm <- lsmeans(tw_rm, ~ Valence * Task)
expect_identical(est_name_from_call(tw_int_lsm), "M")
ow_pairs_lsm <- pairs(ow_me_lsm)
expect_identical(est_name_from_call(ow_pairs_lsm), "\\Delta M")
ow_pairs_lsm_2 <- contrast(ow_me_lsm, interaction = "pairwise")
expect_identical(est_name_from_call(ow_pairs_lsm_2), "\\Delta M")
## Bug reported by shirdekel, #456
tw_pairs_lsm <- pairs(tw_int_lsm)
expect_identical(est_name_from_call(tw_pairs_lsm), "\\Delta M")
tw_pairs_lsm_2 <- contrast(tw_int_lsm, interaction = "pairwise")
expect_identical(est_name_from_call(tw_pairs_lsm_2), "\\Delta M")
tw_int_emm <- emmeans(tw_rm, ~ Valence * Task)
tw_pairs_emm <- pairs(tw_int_emm)
expect_identical(est_name_from_call(tw_pairs_emm), "\\Delta M")
tw_pairs_emm_2 <- contrast(tw_int_emm, interaction = "pairwise")
expect_identical(est_name_from_call(tw_pairs_emm_2), "\\Delta M")
## Univeriate EMM
afex::afex_options(emmeans_model = "univariate")
uni_tw_me_emm <- emmeans(tw_rm, ~ Valence)
expect_identical(est_name_from_call(uni_tw_me_emm), "M")
uni_tw_me_consec_emm <- emmeans(tw_rm, consec ~ Valence)
expect_identical(est_name_from_call(uni_tw_me_consec_emm$emmeans), "M")
expect_identical(est_name_from_call(uni_tw_me_consec_emm$contrasts), "\\Delta M")
uni_tw_pairs_emm <- pairs(emmeans(tw_rm, ~ Valence))
expect_identical(est_name_from_call(uni_tw_pairs_emm), "\\Delta M")
tw_me_emm <- emmeans(tw_rm, ~ Valence* Task)
tw_pairs_contrasts_emm <- contrast(tw_me_emm, "consec", simple = "each")
expect_identical(est_name_from_call(tw_pairs_contrasts_emm$`simple contrasts for Valence`), "\\Delta M")
tw_pairs_contrasts_emm2 <- contrast(tw_me_emm, "consec", simple = "each", combine = TRUE)
expect_identical(est_name_from_call(tw_pairs_contrasts_emm2), "\\Delta M")
## Multivariate EMM
afex::afex_options(emmeans_model = "multivariate")
mw_tw_me_emm <- emmeans(tw_rm, ~ Valence)
expect_identical(est_name_from_call(uni_tw_me_emm), "M")
mw_tw_pairs_emm <- pairs(emmeans(tw_rm, ~ Valence))
expect_identical(est_name_from_call(mw_tw_pairs_emm), "\\Delta M")
# GLM(M)
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
lm_emm <- emmeans(lm.D9, pairwise~group)
expect_identical(est_name_from_call(lm_emm$emmeans), "M")
expect_identical(est_name_from_call(lm_emm$contrasts), "\\Delta M")
glmm <- lme4::glmer(
cbind(incidence, size - incidence) ~ period + (1 | herd)
, data = lme4::cbpp
, family = binomial(link = "logit")
)
glmm_emm_link <- emmeans(glmm, ~ period, type = "link")
expect_identical(est_name_from_call(glmm_emm_link), "\\mathrm{logit}(p)")
glmm_pairs_link <- pairs(glmm_emm_link)
expect_identical(est_name_from_call(glmm_pairs_link), "\\log(\\mathit{OR})")
glmm_emm_resp <- emmeans(glmm, ~ period, type = "response")
expect_identical(est_name_from_call(glmm_emm_resp), "p")
glmm_pairs_resp <- pairs(glmm_emm_resp)
expect_identical(est_name_from_call(glmm_pairs_resp), "\\mathit{OR}")
glmm_pairs_resp2 <- pairs(glmm_emm_resp, ratios = FALSE)
expect_identical(est_name_from_call(glmm_pairs_resp2), "\\Delta \\mathrm{logit}(p)")
glmm2 <- lme4::glmer(
cbind(incidence, size - incidence) ~ period + (1 | herd)
, data = lme4::cbpp
, family = binomial(link = "probit")
)
glmm2_emm_link <- emmeans(glmm2, ~ period, type = "link")
expect_identical(est_name_from_call(glmm2_emm_link), "\\Phi^{-1}(p)")
glmm2_pairs_link <- pairs(glmm2_emm_link)
expect_identical(est_name_from_call(glmm2_pairs_link), "\\Delta \\Phi^{-1}(p)")
glmm2_emm_resp <- emmeans(glmm2, ~ period, type = "response")
expect_identical(est_name_from_call(glmm2_emm_resp), "p")
glmm2_pairs_resp <- pairs(glmm2_emm_resp)
expect_identical(est_name_from_call(glmm2_pairs_resp), "\\Delta \\Phi^{-1}(p)")
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm_mod <- glm(counts ~ outcome + treatment, family = poisson())
glm_emm_link <- emmeans(glm_mod, ~ treatment, type = "link")
expect_identical(est_name_from_call(glm_emm_link), "\\log(M)")
glm_pairs_link <- pairs(glm_emm_link)
expect_identical(est_name_from_call(glm_pairs_link), "\\Delta \\log(M)")
glm_emm_resp <- emmeans(glm_mod, ~ treatment, type = "response")
expect_identical(est_name_from_call(glm_emm_resp), "M")
glm_pairs_resp <- pairs(glm_emm_resp)
expect_identical(est_name_from_call(glm_pairs_resp), "M_{i}/M_{j}")
glm_pairs_resp2 <- pairs(glm_emm_resp, ratio = FALSE)
expect_identical(est_name_from_call(glm_pairs_link), "\\Delta \\log(M)")
}
)
test_that(
"Multiplicity adjustment notes"
, {
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
))
ow_me_emm <- emmeans(tw_rm, ~ Valence | Task)
ow_me_emm_bonf <- summary(ow_me_emm, infer = TRUE, adjust = "bonferroni")
ow_me_emm_bonf_res <- apa_print(ow_me_emm_bonf)
expect_equal(
variable_label(ow_me_emm_bonf_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Bonferroni(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Bonferroni(3)}$"
)
)
ow_me_emm_tukey <- pairs(ow_me_emm, infer = TRUE)
ow_me_emm_tukey_res <- apa_print(ow_me_emm_tukey)
expect_equal(
variable_label(ow_me_emm_tukey_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Tukey(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Tukey(3)}$"
)
)
ow_me_emm_holm <- contrast(ow_me_emm, infer = TRUE, adjust = "holm")
ow_me_emm_holm_res <- apa_print(ow_me_emm_holm)
expect_equal(
variable_label(ow_me_emm_holm_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Bonferroni(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Holm(3)}$"
)
)
ow_me_emm_fdr <- contrast(ow_me_emm, infer = TRUE, adjust = "fdr")
ow_me_emm_fdr_res <- apa_print(ow_me_emm_fdr)
expect_equal(
variable_label(ow_me_emm_fdr_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Bonferroni(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle FDR(3)}$"
)
)
ow_me_emm_dun <- contrast(ow_me_emm, infer = TRUE, adjust = "dunnettx")
ow_me_emm_dun_res <- apa_print(ow_me_emm_dun)
expect_equal(
variable_label(ow_me_emm_dun_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Dunnett(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Dunnett(3)}$"
)
)
ow_me_emm_scheffe <- summary(ow_me_emm, infer = TRUE, adjust = "scheffe")
ow_me_emm_scheffe_res <- apa_print(ow_me_emm_scheffe)
expect_equal(
variable_label(ow_me_emm_scheffe_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
)
)
ow_me_emm_scheffe <- emmeans(ow_me_emm, 1 ~ Valence | Task, adjust = "scheffe")
ow_me_emm_scheffe_res <- apa_print(ow_me_emm_scheffe)
expect_equal(
variable_label(ow_me_emm_scheffe_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
)
)
ow_me_emm_scheffe <- pairs(ow_me_emm, infer = TRUE, adjust = "scheffe")
ow_me_emm_scheffe_res <- apa_print(ow_me_emm_scheffe)
expect_equal(
variable_label(ow_me_emm_scheffe_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Scheff\\'e(2)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Scheff\\'e(2)}$"
)
)
family_mark_reference <- letters[
rep(1:nlevels(tw_rm_data$Task), each = nlevels(tw_rm_data$Valence))
]
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$conf.int
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$adj.p.value
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
ow_me_emm <- emmeans(tw_rm, ~ Task | Valence)
ow_me_emm_bonf <- summary(ow_me_emm, infer = TRUE, adjust = "bonferroni")
ow_me_emm_bonf_res <- apa_print(ow_me_emm_bonf)
family_mark_reference <- letters[
rep(1:nlevels(tw_rm_data$Valence), each = nlevels(tw_rm_data$Task))
]
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$conf.int
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$adj.p.value
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
df <- data.frame(
errors = floor(runif(n=320,min=0,max=30))
, session=c("t1","t2","t3","t4")
)
glmm <- glm(errors ~ session, df, family = "poisson")
glmm_pairs <- emmeans(glmm, pairwise~session, type="response")
glmm_pairs_res <- apa_print(glmm_pairs$contrasts)
expect_false(any(grepl("[a-z]", glmm_pairs_res$table$conf.int)))
expect_false(any(grepl("[a-z]", glmm_pairs_res$table$p.value)))
}
)
test_that(
"Regression"
, {
# Typesetting of numeric predictor values
# https://github.com/crsh/papaja/issues/445
iris$Sepal.Length <- iris$Sepal.Length * 1000
my_lm <- lm(
Sepal.Width ~ Sepal.Length + Petal.Width + Petal.Length
, data = iris
)
ats <- c(1000, 2000, 3000)
my_lm_emm <- emmeans::emmeans(
my_lm
, ~Sepal.Length
, at = list(Sepal.Length = ats)
)
my_lm_emm_output <- apa_print(my_lm_emm)
expect_apa_results(
my_lm_emm_output
, labels = list(
Sepal.Length = "Sepal.Length"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = c("X1000", "X2000", "X3000")
)
expect_apa_term(
my_lm_emm_output
, term = "X1000"
, estimate = "$M = 0.12$, 95\\% CI $[-0.48, 0.71]$"
, statistic = "$t(146) = 0.39$, $p = .699$"
)
expect_true(
all(my_lm_emm_output$table$Sepal.Length == apa_num(ats))
)
# https://github.com/crsh/papaja/issues/200
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
emm_glm <- emmeans::emmeans(glm.D93, pairwise~treatment)
emm_glm_output <- apa_print(emm_glm$emmeans)
expect_apa_results(
emm_glm_output
, labels = list(
treatment = "treatment"
, estimate = "$\\log(M)$"
, conf.int = "95\\% CI"
, statistic = "$z$"
, p.value = "$p$"
)
, term_names = c("X1", "X2", "X3")
)
expect_apa_term(
emm_glm_output
, term = "X1"
, estimate = "$\\log(M) = 2.80$, 95\\% CI $[2.52, 3.07]$"
, statistic = "$z = 19.65$, $p < .001$"
)
# pigs_lm <- lm(log(conc) ~ source * percent, data = emmeans::pigs)
#
# pigs_lm_emm <- emmeans::emmeans(pigs_lm, ~ source)
# pigs_pairs_emm_output <- apa_print(pairs(pigs_lm_emm, type = "response"))
#
# # table --------------------------------------------------------------------
# expect_identical(
# object = pigs_pairs_emm_output$table$estimate
# , expected = structure(
# c("0.77", "0.66", "0.86")
# , label = "Ratio"
# , class = c("tiny_labelled", "character")
# )
# )
#
#
# noise.lm <- lm(noise ~ size * type * side, data = emmeans::auto.noise)
# noise.emm <- emmeans::emmeans(noise.lm, ~ size * side * type)
# apa_print(emmeans::contrast(noise.emm, "consec", simple = "each", combine = TRUE, adjust = "mvt")) # Table order is fucked up
#
#
# org.int <- lm(cbind(sales1, sales2) ~ price1 * price2 + day + store, data = emmeans::oranges)
# emmeans::emtrends(org.int, ~ variety, var = "price1", mult.name = "variety")
# emmeans::emtrends(org.int, pairwise ~ variety, var = "price1", mult.name = "variety")
#
}
)
test_that(
"emtrends"
, {
skip("emtrends() is not yet supported.")
library("emmeans")
cars_lm <- lm(mpg ~ qsec, data = mtcars)
summary(cars_lm)
cars_em <- emtrends(cars_lm, ~ 1, var = "qsec")
onesided <- summary(cars_em, infer = TRUE, side = ">")
apa_print(onesided)
# broom::tidy(summary(cars_em, infer = TRUE, side = ">"))
}
)
crsh/papaja documentation built on Feb. 21, 2024, 6:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
context("apa_print() for emmeans/lsmeans")
test_that(
"ANOVA"
, {
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
))
# Main effect ------------------------------------------------------
tw_me_emm <- emmeans::emmeans(tw_rm, ~ Valence)
tw_me_emm_output <- expect_warning(
apa_print(tw_me_emm, conf.l = .95)
, regexp = "Using argument 'conf.level' in calls to 'apa_print()' is deprecated. Please use 'conf.int' instead."
, fixed = TRUE
)
expect_apa_results(
tw_me_emm_output
, labels = list(
Valence = "Valence"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = levels(tw_rm_data$Valence)
)
expect_apa_term(
tw_me_emm_output
, term = "Neg"
, estimate = "$M = 11.00$, 95\\% CI $[6.34, 15.66]$"
, statistic = "$t(4.63) = 6.21$, $p = .002$"
)
# Alternative calls
tw_me_lsm <- emmeans::lsmeans(tw_rm, ~ Valence)
tw_me_lsm_output <- apa_print(tw_me_lsm)
tw_me_lsm_output2 <- apa_print(
summary(tw_me_lsm, infer = TRUE)
, est_name = "M"
)
tw_me_emm_output2 <- apa_print(
summary(tw_me_emm, infer = TRUE)
, est_name = "M"
)
expect_identical(tw_me_emm_output, tw_me_emm_output2)
expect_identical(tw_me_lsm_output, tw_me_lsm_output2)
expect_identical(tw_me_lsm_output, tw_me_emm_output2)
expect_warning(
tw_me_emm_p_output <- apa_print(
summary(tw_me_emm, infer = c(FALSE, TRUE))
, est_name = "M"
)
)
expect_identical(
tw_me_lsm_output$table[, c("Valence", "estimate", "statistic", "df", "p.value")]
, tw_me_emm_p_output$table
)
tw_me_emm_ci_output <- apa_print(
summary(tw_me_emm, infer = c(TRUE, FALSE))
, est_name = "M"
)
expect_apa_term(
tw_me_emm_ci_output
, term = "Neg"
, estimate = "$M = 11.00$, 95\\% CI $[6.34, 15.66]$"
, statistic = NULL
)
expect_identical(
tw_me_lsm_output$table[, c("Valence", "estimate", "conf.int")]
, tw_me_emm_ci_output$table
)
expect_apa_results(
apa_print(summary(tw_me_lsm, infer = c(FALSE, FALSE)))
, labels = list(
Valence = "Valence"
, estimate = "$\\hat{\\theta}$"
)
, term_names = levels(tw_rm_data$Valence)
)
tw_me_emm_output3 <- apa_print(
summary(tw_me_emm, infer = TRUE, side = ">")
)
# https://github.com/crsh/papaja/issues/456#issuecomment-901653372
expect_apa_term(
tw_me_emm_output3
, term = "Neg"
, estimate = "$\\hat{\\theta} = 11.00$, 95\\% CI $[7.37, \\infty]$"
, statistic = "$t(4.63) = 6.21$, $p = .001$"
)
# Interaction ------------------------------------------------------
tw_int_emm <- emmeans::emmeans(tw_rm, ~ Valence * Task)
tw_int_emm_output <- apa_print(tw_int_emm)
term_names <- apply(expand.grid(levels(tw_rm_data$Valence), levels(tw_rm_data$Task)), 1, paste, collapse = "_")
expect_apa_results(
tw_int_emm_output
, labels = list(
Task = "Task"
, Valence = "Valence"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = term_names
)
expect_apa_term(
tw_int_emm_output
, term = "Neg_Cued"
, estimate = "$M = 11.80$, 95\\% CI $[7.17, 16.43]$"
, statistic = "$t(5.52) = 6.37$, $p < .001$"
)
# Alternative calls
tw_int_lsm <- emmeans::lsmeans(tw_rm, ~ Valence * Task)
tw_int_lsm_output <- apa_print(tw_int_lsm)
tw_int_lsm_output2 <- apa_print(
summary(tw_int_lsm, infer = TRUE)
, est_name = "M"
)
tw_int_emm_output2 <- apa_print(
summary(tw_int_emm, infer = TRUE)
, est_name = "M"
)
expect_identical(tw_int_emm_output, tw_int_emm_output2)
expect_identical(tw_int_lsm_output, tw_int_lsm_output2)
expect_identical(tw_int_lsm_output, tw_int_emm_output2)
# Simple effects ---------------------------------------------------
tw_se_emm <- emmeans::emmeans(tw_rm, ~ Valence | Task)
tw_se_emm_output <- apa_print(tw_se_emm)
expect_apa_results(
tw_se_emm_output
, labels = list(
Task = "Task"
, Valence = "Valence"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = term_names
)
# Sort
tw_int_emm_output$estimate <- tw_int_emm_output$estimate[names(tw_se_emm_output$estimate)]
tw_int_emm_output$statistic <- tw_int_emm_output$statistic[names(tw_int_emm_output$estimate)]
tw_int_emm_output$full_result <- tw_int_emm_output$full_result[names(tw_int_emm_output$estimate)]
tw_int_emm_output$table <- tw_int_emm_output$table[rownames(tw_se_emm_output$table), ]
expect_identical(tw_se_emm_output, tw_int_emm_output)
# Complex output (four factors) ------------------------------------
data(obk.long, package = "afex")
fw_mixed <- suppressWarnings(afex::aov_ez(
id = "id"
, dv = "value"
, data = obk.long
, between = c("treatment", "gender")
, within = c("phase", "hour")
, observed = "gender"
))
fw_mixed_emm <- emmeans::emmeans(fw_mixed, ~ phase * hour | treatment * gender)
fw_mixed_emm_output <- apa_print(fw_mixed_emm)
expect_apa_results(
fw_mixed_emm_output
, labels = list(
gender = "gender"
, treatment = "treatment"
, hour = "hour"
, phase = "phase"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
)
# table ------------------------------------------------------------
expect_identical(
object = fw_mixed_emm_output$table$gender
, expected = structure(
c(
"F"
, rep("", nrow(fw_mixed_emm_output$table) / 2 - 1)
, "M"
, rep("", nrow(fw_mixed_emm_output$table) / 2 - 1)
)
, class = c("tiny_labelled", "character")
, label = "gender"
)
)
expect_identical(
object = fw_mixed_emm_output$table$treatment
, expected = structure(
rep(
c(
"Control"
, rep("", nrow(fw_mixed_emm_output$table) / 6 - 1)
, "A"
, rep("", nrow(fw_mixed_emm_output$table) / 6 - 1)
, "B"
, rep("", nrow(fw_mixed_emm_output$table) / 6 - 1)
)
, 2
)
, class = c("tiny_labelled", "character")
, label = "treatment"
)
)
expect_identical(
object = fw_mixed_emm_output$table$hour
, expected = structure(
rep(
c(
"X1"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X2"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X3"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X4"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
, "X5"
, rep("", nrow(fw_mixed_emm_output$table) / 30 - 1)
)
, 6
)
, class = c("tiny_labelled", "character")
, label = "hour"
)
)
expect_identical(
object = fw_mixed_emm_output$table$phase
, expected = structure(
rep(
c("Fup", "Post", "Pre")
, 30
)
, class = c("tiny_labelled", "character")
, label = "phase"
)
)
# Alternative calls
fw_mixed_lsm <- emmeans::lsmeans(fw_mixed, ~ phase * hour | treatment * gender)
fw_mixed_lsm_output <- apa_print(fw_mixed_lsm)
fw_mixed_lsm_output2 <- apa_print(
summary(fw_mixed_lsm, infer = TRUE)
, est_name = "M"
)
fw_mixed_emm_output2 <- apa_print(
summary(fw_mixed_emm, infer = TRUE)
, est_name = "M"
)
expect_identical(fw_mixed_lsm_output, fw_mixed_lsm_output2)
expect_identical(fw_mixed_emm_output, fw_mixed_emm_output2)
expect_identical(fw_mixed_lsm_output, fw_mixed_emm_output2)
# Pairs ------------------------------------------------------------
tw_me_pairs_emm <- pairs(tw_me_emm)
tw_me_pairs_emm_output <- apa_print(tw_me_pairs_emm)
expect_apa_results(
tw_me_pairs_emm_output
, labels = list(
contrast = "Contrast"
, estimate = "$\\Delta M$"
, conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Tukey(3)}$"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Tukey(3)}$"
)
, term_names = c("Neg_Neu", "Neg_Pos", "Neu_Pos")
)
expect_apa_term(
tw_me_pairs_emm_output
, term = "Neg_Neu"
, estimate = "$\\Delta M = -1.10$, $95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Tukey(3)}$ $[-3.44, 1.24]$"
, statistic = "$t(8) = -1.34$, $p_\\mathrm{\\scriptstyle Tukey(3)} = .413$"
)
## Custom contrast names
tw_me_pairs_emm_output <- apa_print(
tw_me_pairs_emm
, contrast_names = letters[1:3]
)
expect_identical(
object = tw_me_pairs_emm_output$table$contrast
, expected = structure(
letters[1:3]
, label = "Contrast"
, class = c("tiny_labelled", "character")
)
)
# Custom contrasts
tw_se_contrast_emm <- emmeans::contrast(
tw_se_emm
, method = list("Positive - Negative" = c(1, 0, -1))
)
tw_se_contrast_emm_output <- apa_print(tw_se_contrast_emm)
tw_se_contrast_emm_output2 <- apa_print(
summary(tw_se_contrast_emm, infer = TRUE)
, est_name = "\\Delta M"
)
expect_identical(
tw_se_contrast_emm_output
, tw_se_contrast_emm_output2
)
# Simple contrasts
simple_pairs <- pairs(
emmeans::emmeans(fw_mixed, ~ gender | treatment)
)
tw_between_emm <- emmeans::emmeans(fw_mixed, ~ treatment * gender)
simple_pairs2 <- pairs(
tw_between_emm
, simple = "gender"
)
simple_pairs_output <- apa_print(simple_pairs)
simple_pairs2_output <- apa_print(simple_pairs2)
expect_apa_results(
simple_pairs_output
, term_names = c("Control_F_M", "A_F_M", "B_F_M")
)
expect_identical(simple_pairs_output, simple_pairs2_output)
simple_contrasts <- emmeans::contrast(
tw_between_emm
, "consec"
, simple = "each"
, combine = TRUE
, adjust = "none"
)
simple_contrasts_output <- apa_print(simple_contrasts)
simple_contrasts2 <- emmeans::contrast(
tw_between_emm
, "consec"
, simple = "each"
, adjust = "none"
)
simple_contrasts2_output1 <- apa_print(simple_contrasts2$`simple contrasts for gender`)
simple_contrasts2_output2 <- apa_print(simple_contrasts2$`simple contrasts for treatment`)
expect_identical(
simple_contrasts_output$full_result[grepl("M_F$", names(simple_contrasts_output$full_result))]
, simple_contrasts2_output1$full_result
)
expect_identical(
simple_contrasts_output$full_result[grepl("^(F|M)_", names(simple_contrasts_output$full_result))]
, simple_contrasts2_output2$full_result
)
}
)
test_that(
"Joint tests"
, {
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
# Joint tests
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
, anova_table = list(correction = "none")
))
tw_rm_output <- suppressWarnings(apa_print(tw_rm))
tw_rm_emm <- emmeans::emmeans(tw_rm$aov, ~ Task * Valence)
## All terms
emm_aov <- emmeans::joint_tests(tw_rm_emm)
emm_aov_output <- apa_print(emm_aov)
expect_apa_results(
emm_aov_output
, labels = list(
term = "Effect"
, statistic = "$F$"
, df = "$\\mathit{df}$"
, df.residual = "$\\mathit{df}_{\\mathrm{res}}$"
, p.value = "$p$"
)
, term_names = names(tw_rm_output$estimate)
)
expect_apa_term(
emm_aov_output
, term = "Task"
, estimate = NULL
, statistic = tw_rm_output$statistic$Task
)
## Split by
emm_split_aov <- emmeans::joint_tests(tw_rm_emm, by = "Task")
emm_split_aov_output <- apa_print(emm_split_aov)
expect_apa_results(
emm_split_aov_output
, labels = list(
Task = "Task"
, term = "Effect"
, statistic = "$F$"
, df = "$\\mathit{df}$"
, df.residual = "$\\mathit{df}_{\\mathrm{res}}$"
, p.value = "$p$"
)
, term_names = paste("Valence", levels(tw_rm_data$Task), sep = "_")
, table_terms = rep("Valence", 2) # Levels split by are in "Task" column
)
expect_apa_term(
emm_split_aov_output
, term = "Valence_Cued"
, estimate = NULL
, statistic = "$F(2, 15.58) = 1.46$, $p = .263$"
)
# Ensure proper sorting of terms
load("data/mixed_data.rdata")
unsorted_aov <- afex::aov_4(formula = Recall ~ Gender * Dosage * (Task * Valence |Subject), data = mixed_data, fun_aggregate = mean)
unsorted_emm <- emmeans::joint_tests(unsorted_aov, by = "Gender")
apa_out <- apa_print(unsorted_emm)
expect_apa_results(
apa_out
, labels = list(
Gender = "Gender"
, term = "Effect"
, statistic = "$F$"
, df = "$\\mathit{df}$"
, df.residual = "$\\mathit{df}_{\\mathrm{res}}$"
, p.value = "$p$"
)
, term_names = papaja:::sanitize_terms(paste(unlabel(gsub(apa_out$table$term, pattern = " $\\times$ ", replacement = "_", fixed = TRUE)), apa_out$table$Gender, sep = "_"))
, table_terms = beautify_terms(data.frame(unsorted_emm)$model.term)
)
}
)
test_that(
"Estimate name guessing"
, {
# ANOVA
library("emmeans")
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
))
ow_me_lsm <- lsmeans(tw_rm, ~ Valence)
expect_identical(est_name_from_call(ow_me_lsm), "M")
tw_int_lsm <- lsmeans(tw_rm, ~ Valence * Task)
expect_identical(est_name_from_call(tw_int_lsm), "M")
ow_pairs_lsm <- pairs(ow_me_lsm)
expect_identical(est_name_from_call(ow_pairs_lsm), "\\Delta M")
ow_pairs_lsm_2 <- contrast(ow_me_lsm, interaction = "pairwise")
expect_identical(est_name_from_call(ow_pairs_lsm_2), "\\Delta M")
## Bug reported by shirdekel, #456
tw_pairs_lsm <- pairs(tw_int_lsm)
expect_identical(est_name_from_call(tw_pairs_lsm), "\\Delta M")
tw_pairs_lsm_2 <- contrast(tw_int_lsm, interaction = "pairwise")
expect_identical(est_name_from_call(tw_pairs_lsm_2), "\\Delta M")
tw_int_emm <- emmeans(tw_rm, ~ Valence * Task)
tw_pairs_emm <- pairs(tw_int_emm)
expect_identical(est_name_from_call(tw_pairs_emm), "\\Delta M")
tw_pairs_emm_2 <- contrast(tw_int_emm, interaction = "pairwise")
expect_identical(est_name_from_call(tw_pairs_emm_2), "\\Delta M")
## Univeriate EMM
afex::afex_options(emmeans_model = "univariate")
uni_tw_me_emm <- emmeans(tw_rm, ~ Valence)
expect_identical(est_name_from_call(uni_tw_me_emm), "M")
uni_tw_me_consec_emm <- emmeans(tw_rm, consec ~ Valence)
expect_identical(est_name_from_call(uni_tw_me_consec_emm$emmeans), "M")
expect_identical(est_name_from_call(uni_tw_me_consec_emm$contrasts), "\\Delta M")
uni_tw_pairs_emm <- pairs(emmeans(tw_rm, ~ Valence))
expect_identical(est_name_from_call(uni_tw_pairs_emm), "\\Delta M")
tw_me_emm <- emmeans(tw_rm, ~ Valence* Task)
tw_pairs_contrasts_emm <- contrast(tw_me_emm, "consec", simple = "each")
expect_identical(est_name_from_call(tw_pairs_contrasts_emm$`simple contrasts for Valence`), "\\Delta M")
tw_pairs_contrasts_emm2 <- contrast(tw_me_emm, "consec", simple = "each", combine = TRUE)
expect_identical(est_name_from_call(tw_pairs_contrasts_emm2), "\\Delta M")
## Multivariate EMM
afex::afex_options(emmeans_model = "multivariate")
mw_tw_me_emm <- emmeans(tw_rm, ~ Valence)
expect_identical(est_name_from_call(uni_tw_me_emm), "M")
mw_tw_pairs_emm <- pairs(emmeans(tw_rm, ~ Valence))
expect_identical(est_name_from_call(mw_tw_pairs_emm), "\\Delta M")
# GLM(M)
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
lm_emm <- emmeans(lm.D9, pairwise~group)
expect_identical(est_name_from_call(lm_emm$emmeans), "M")
expect_identical(est_name_from_call(lm_emm$contrasts), "\\Delta M")
glmm <- lme4::glmer(
cbind(incidence, size - incidence) ~ period + (1 | herd)
, data = lme4::cbpp
, family = binomial(link = "logit")
)
glmm_emm_link <- emmeans(glmm, ~ period, type = "link")
expect_identical(est_name_from_call(glmm_emm_link), "\\mathrm{logit}(p)")
glmm_pairs_link <- pairs(glmm_emm_link)
expect_identical(est_name_from_call(glmm_pairs_link), "\\log(\\mathit{OR})")
glmm_emm_resp <- emmeans(glmm, ~ period, type = "response")
expect_identical(est_name_from_call(glmm_emm_resp), "p")
glmm_pairs_resp <- pairs(glmm_emm_resp)
expect_identical(est_name_from_call(glmm_pairs_resp), "\\mathit{OR}")
glmm_pairs_resp2 <- pairs(glmm_emm_resp, ratios = FALSE)
expect_identical(est_name_from_call(glmm_pairs_resp2), "\\Delta \\mathrm{logit}(p)")
glmm2 <- lme4::glmer(
cbind(incidence, size - incidence) ~ period + (1 | herd)
, data = lme4::cbpp
, family = binomial(link = "probit")
)
glmm2_emm_link <- emmeans(glmm2, ~ period, type = "link")
expect_identical(est_name_from_call(glmm2_emm_link), "\\Phi^{-1}(p)")
glmm2_pairs_link <- pairs(glmm2_emm_link)
expect_identical(est_name_from_call(glmm2_pairs_link), "\\Delta \\Phi^{-1}(p)")
glmm2_emm_resp <- emmeans(glmm2, ~ period, type = "response")
expect_identical(est_name_from_call(glmm2_emm_resp), "p")
glmm2_pairs_resp <- pairs(glmm2_emm_resp)
expect_identical(est_name_from_call(glmm2_pairs_resp), "\\Delta \\Phi^{-1}(p)")
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm_mod <- glm(counts ~ outcome + treatment, family = poisson())
glm_emm_link <- emmeans(glm_mod, ~ treatment, type = "link")
expect_identical(est_name_from_call(glm_emm_link), "\\log(M)")
glm_pairs_link <- pairs(glm_emm_link)
expect_identical(est_name_from_call(glm_pairs_link), "\\Delta \\log(M)")
glm_emm_resp <- emmeans(glm_mod, ~ treatment, type = "response")
expect_identical(est_name_from_call(glm_emm_resp), "M")
glm_pairs_resp <- pairs(glm_emm_resp)
expect_identical(est_name_from_call(glm_pairs_resp), "M_{i}/M_{j}")
glm_pairs_resp2 <- pairs(glm_emm_resp, ratio = FALSE)
expect_identical(est_name_from_call(glm_pairs_link), "\\Delta \\log(M)")
}
)
test_that(
"Multiplicity adjustment notes"
, {
emm_basis.afex_aov <- afex:::emm_basis.afex_aov
load("data/tw_rm_data.rdata")
tw_rm <- suppressWarnings(afex::aov_ez(
data = tw_rm_data
, id = "Subject"
, dv = "Recall"
, within = c("Task", "Valence")
))
ow_me_emm <- emmeans(tw_rm, ~ Valence | Task)
ow_me_emm_bonf <- summary(ow_me_emm, infer = TRUE, adjust = "bonferroni")
ow_me_emm_bonf_res <- apa_print(ow_me_emm_bonf)
expect_equal(
variable_label(ow_me_emm_bonf_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Bonferroni(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Bonferroni(3)}$"
)
)
ow_me_emm_tukey <- pairs(ow_me_emm, infer = TRUE)
ow_me_emm_tukey_res <- apa_print(ow_me_emm_tukey)
expect_equal(
variable_label(ow_me_emm_tukey_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Tukey(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Tukey(3)}$"
)
)
ow_me_emm_holm <- contrast(ow_me_emm, infer = TRUE, adjust = "holm")
ow_me_emm_holm_res <- apa_print(ow_me_emm_holm)
expect_equal(
variable_label(ow_me_emm_holm_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Bonferroni(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Holm(3)}$"
)
)
ow_me_emm_fdr <- contrast(ow_me_emm, infer = TRUE, adjust = "fdr")
ow_me_emm_fdr_res <- apa_print(ow_me_emm_fdr)
expect_equal(
variable_label(ow_me_emm_fdr_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Bonferroni(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle FDR(3)}$"
)
)
ow_me_emm_dun <- contrast(ow_me_emm, infer = TRUE, adjust = "dunnettx")
ow_me_emm_dun_res <- apa_print(ow_me_emm_dun)
expect_equal(
variable_label(ow_me_emm_dun_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Dunnett(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Dunnett(3)}$"
)
)
ow_me_emm_scheffe <- summary(ow_me_emm, infer = TRUE, adjust = "scheffe")
ow_me_emm_scheffe_res <- apa_print(ow_me_emm_scheffe)
expect_equal(
variable_label(ow_me_emm_scheffe_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
)
)
ow_me_emm_scheffe <- emmeans(ow_me_emm, 1 ~ Valence | Task, adjust = "scheffe")
ow_me_emm_scheffe_res <- apa_print(ow_me_emm_scheffe)
expect_equal(
variable_label(ow_me_emm_scheffe_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Scheff\\'e(3)}$"
)
)
ow_me_emm_scheffe <- pairs(ow_me_emm, infer = TRUE, adjust = "scheffe")
ow_me_emm_scheffe_res <- apa_print(ow_me_emm_scheffe)
expect_equal(
variable_label(ow_me_emm_scheffe_res$table[, c("conf.int", "adj.p.value")])
, list(
conf.int = "$95\\%\\ \\mathrm{CI}_\\mathrm{\\scriptstyle Scheff\\'e(2)}$"
, adj.p.value = "$p_\\mathrm{\\scriptstyle Scheff\\'e(2)}$"
)
)
family_mark_reference <- letters[
rep(1:nlevels(tw_rm_data$Task), each = nlevels(tw_rm_data$Valence))
]
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$conf.int
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$adj.p.value
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
ow_me_emm <- emmeans(tw_rm, ~ Task | Valence)
ow_me_emm_bonf <- summary(ow_me_emm, infer = TRUE, adjust = "bonferroni")
ow_me_emm_bonf_res <- apa_print(ow_me_emm_bonf)
family_mark_reference <- letters[
rep(1:nlevels(tw_rm_data$Valence), each = nlevels(tw_rm_data$Task))
]
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$conf.int
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
family_marks <- .str_extract_first(
.str_extract_first(
ow_me_emm_bonf_res$table$adj.p.value
, "[a-z]\\$$"
)
, "[a-z]"
)
family_marks <- unlabel(family_marks)
expect_equal(
family_marks
, family_mark_reference
)
df <- data.frame(
errors = floor(runif(n=320,min=0,max=30))
, session=c("t1","t2","t3","t4")
)
glmm <- glm(errors ~ session, df, family = "poisson")
glmm_pairs <- emmeans(glmm, pairwise~session, type="response")
glmm_pairs_res <- apa_print(glmm_pairs$contrasts)
expect_false(any(grepl("[a-z]", glmm_pairs_res$table$conf.int)))
expect_false(any(grepl("[a-z]", glmm_pairs_res$table$p.value)))
}
)
test_that(
"Regression"
, {
# Typesetting of numeric predictor values
# https://github.com/crsh/papaja/issues/445
iris$Sepal.Length <- iris$Sepal.Length * 1000
my_lm <- lm(
Sepal.Width ~ Sepal.Length + Petal.Width + Petal.Length
, data = iris
)
ats <- c(1000, 2000, 3000)
my_lm_emm <- emmeans::emmeans(
my_lm
, ~Sepal.Length
, at = list(Sepal.Length = ats)
)
my_lm_emm_output <- apa_print(my_lm_emm)
expect_apa_results(
my_lm_emm_output
, labels = list(
Sepal.Length = "Sepal.Length"
, estimate = "$M$"
, conf.int = "95\\% CI"
, statistic = "$t$"
, df = "$\\mathit{df}$"
, p.value = "$p$"
)
, term_names = c("X1000", "X2000", "X3000")
)
expect_apa_term(
my_lm_emm_output
, term = "X1000"
, estimate = "$M = 0.12$, 95\\% CI $[-0.48, 0.71]$"
, statistic = "$t(146) = 0.39$, $p = .699$"
)
expect_true(
all(my_lm_emm_output$table$Sepal.Length == apa_num(ats))
)
# https://github.com/crsh/papaja/issues/200
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
emm_glm <- emmeans::emmeans(glm.D93, pairwise~treatment)
emm_glm_output <- apa_print(emm_glm$emmeans)
expect_apa_results(
emm_glm_output
, labels = list(
treatment = "treatment"
, estimate = "$\\log(M)$"
, conf.int = "95\\% CI"
, statistic = "$z$"
, p.value = "$p$"
)
, term_names = c("X1", "X2", "X3")
)
expect_apa_term(
emm_glm_output
, term = "X1"
, estimate = "$\\log(M) = 2.80$, 95\\% CI $[2.52, 3.07]$"
, statistic = "$z = 19.65$, $p < .001$"
)
# pigs_lm <- lm(log(conc) ~ source * percent, data = emmeans::pigs)
#
# pigs_lm_emm <- emmeans::emmeans(pigs_lm, ~ source)
# pigs_pairs_emm_output <- apa_print(pairs(pigs_lm_emm, type = "response"))
#
# # table --------------------------------------------------------------------
# expect_identical(
# object = pigs_pairs_emm_output$table$estimate
# , expected = structure(
# c("0.77", "0.66", "0.86")
# , label = "Ratio"
# , class = c("tiny_labelled", "character")
# )
# )
#
#
# noise.lm <- lm(noise ~ size * type * side, data = emmeans::auto.noise)
# noise.emm <- emmeans::emmeans(noise.lm, ~ size * side * type)
# apa_print(emmeans::contrast(noise.emm, "consec", simple = "each", combine = TRUE, adjust = "mvt")) # Table order is fucked up
#
#
# org.int <- lm(cbind(sales1, sales2) ~ price1 * price2 + day + store, data = emmeans::oranges)
# emmeans::emtrends(org.int, ~ variety, var = "price1", mult.name = "variety")
# emmeans::emtrends(org.int, pairwise ~ variety, var = "price1", mult.name = "variety")
#
}
)
test_that(
"emtrends"
, {
skip("emtrends() is not yet supported.")
library("emmeans")
cars_lm <- lm(mpg ~ qsec, data = mtcars)
summary(cars_lm)
cars_em <- emtrends(cars_lm, ~ 1, var = "qsec")
onesided <- summary(cars_em, infer = TRUE, side = ">")
apa_print(onesided)
# broom::tidy(summary(cars_em, infer = TRUE, side = ">"))
}
)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.