A3_Q2
library(pacman)
p_load(afex, apa, effectsize, rstatix, qqplotr, car, pwr, rio, gt, psych, gtsummary, tidyverse, ggplot2, apaTables, janitor, ez, performance, emmeans, viridis, report)
dat2 <- import("A3_Q2.sav")
glimpse(dat2)
## Rows: 36
## Columns: 3
## $ noise <dbl> 810, 820, 820, 840, 840, 845, 785, 790, 785, 835, 835, 835, 84…
## $ carsize <dbl> 1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1,…
## $ type <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2,…
dat2 <- dat2 |>
mutate(id = row_number(),
type = as.factor(type),
carsize = as.factor(carsize)) |>
relocate(id, .before = 1)
dat2
## id noise carsize type
## 1 1 810 1 1
## 2 2 820 1 1
## 3 3 820 1 1
## 4 4 840 2 1
## 5 5 840 2 1
## 6 6 845 2 1
## 7 7 785 3 1
## 8 8 790 3 1
## 9 9 785 3 1
## 10 10 835 1 1
## 11 11 835 1 1
## 12 12 835 1 1
## 13 13 845 2 1
## 14 14 855 2 1
## 15 15 850 2 1
## 16 16 760 3 1
## 17 17 760 3 1
## 18 18 770 3 1
## 19 19 820 1 2
## 20 20 820 1 2
## 21 21 820 1 2
## 22 22 820 2 2
## 23 23 820 2 2
## 24 24 825 2 2
## 25 25 775 3 2
## 26 26 775 3 2
## 27 27 775 3 2
## 28 28 825 1 2
## 29 29 825 1 2
## 30 30 825 1 2
## 31 31 815 2 2
## 32 32 825 2 2
## 33 33 825 2 2
## 34 34 770 3 2
## 35 35 760 3 2
## 36 36 765 3 2
Descriptives
describeBy(dat2$noise,
group = dat2$carsize:dat2$type)
##
## Descriptive statistics by group
## group: 1:1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 6 825.83 10.68 827.5 825.83 11.12 810 835 25 -0.28 -1.9 4.36
## ------------------------------------------------------------
## group: 1:2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 6 822.5 2.74 822.5 822.5 3.71 820 825 5 0 -2.31 1.12
## ------------------------------------------------------------
## group: 2:1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 6 845.83 5.85 845 845.83 7.41 840 855 15 0.37 -1.62 2.39
## ------------------------------------------------------------
## group: 2:2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 6 821.67 4.08 822.5 821.67 3.71 815 825 10 -0.48 -1.58 1.67
## ------------------------------------------------------------
## group: 3:1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 6 775 13.42 777.5 775 14.83 760 790 30 -0.1 -2.11 5.48
## ------------------------------------------------------------
## group: 3:2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 6 770 6.32 772.5 770 3.71 760 775 15 -0.49 -1.7 2.58
Visualizing Differences
dat2 |>
ggplot(aes(x = carsize, y = noise)) +
geom_boxplot(aes(color = type)) +
geom_jitter(aes(color = type), width = .25) +
scale_color_viridis(discrete = TRUE) +
theme_bw()
dat2 |>
group_by(carsize, type) |>
summarize(mean = mean(noise)) |>
ggplot(aes(x = carsize, y = mean, color = type, group = type)) +
geom_line()
## `summarise()` has grouped output by 'carsize'. You can override using the
## `.groups` argument.
Model
model2 <- aov_ez(id = "id",
data = dat2,
between = c("carsize", "type"),
dv = "noise",
type = 3,
anova_table = list(es = "pes"))
## Contrasts set to contr.sum for the following variables: carsize, type
check_model(model2)
check_homogeneity(model2)
## Warning: Variances differ between groups (Levene's Test, p = 0.000).
ANOVA
model2_tbl <- model2$lm
apa.aov.table(model2_tbl, conf.level = .95, filename = "assignment_2_anova_tbl.doc")
##
##
## ANOVA results using dv as the dependent variable
##
##
## Predictor SS df MS F p partial_eta2
## (Intercept) 23627700.69 1 23627700.69 361187.78 .000
## carsize 26051.39 2 13025.69 199.12 .000 .93
## type 1056.25 1 1056.25 16.15 .000 .35
## carsize x type 804.17 2 402.08 6.15 .006 .29
## Error 1962.50 30 65.42
## CI_95_partial_eta2
##
## [.86, .95]
## [.09, .55]
## [.03, .48]
##
##
## Note: Values in square brackets indicate the bounds of the 95% confidence interval for partial eta-squared
emmeans(model2, specs = ~carsize*type) |>
contrast(interaction = "pairwise",
simple = "each",
combine = FALSE,
adjust = "tukey")
## $`simple contrasts for carsize`
## type = 1:
## carsize_pairwise estimate SE df t.ratio p.value
## 1 - 2 -20.000 4.67 30 -4.283 0.0005
## 1 - 3 50.833 4.67 30 10.886 <.0001
## 2 - 3 70.833 4.67 30 15.169 <.0001
##
## type = 2:
## carsize_pairwise estimate SE df t.ratio p.value
## 1 - 2 0.833 4.67 30 0.178 0.9826
## 1 - 3 52.500 4.67 30 11.243 <.0001
## 2 - 3 51.667 4.67 30 11.064 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## $`simple contrasts for type`
## carsize = 1:
## type_pairwise estimate SE df t.ratio p.value
## 1 - 2 3.33 4.67 30 0.714 0.4808
##
## carsize = 2:
## type_pairwise estimate SE df t.ratio p.value
## 1 - 2 24.17 4.67 30 5.175 <.0001
##
## carsize = 3:
## type_pairwise estimate SE df t.ratio p.value
## 1 - 2 5.00 4.67 30 1.071 0.2928
library(report)
aov(noise ~ carsize*type, data = dat2) %>%
report()
## The ANOVA (formula: noise ~ carsize * type) suggests that:
##
## - The main effect of carsize is statistically significant and large (F(2, 30) =
## 199.12, p < .001; Eta2 (partial) = 0.93, 95% CI [0.89, 1.00])
## - The main effect of type is statistically significant and large (F(1, 30) =
## 16.15, p < .001; Eta2 (partial) = 0.35, 95% CI [0.13, 1.00])
## - The interaction between carsize and type is statistically significant and
## large (F(2, 30) = 6.15, p = 0.006; Eta2 (partial) = 0.29, 95% CI [0.06, 1.00])
##
## Effect sizes were labelled following Field's (2013) recommendations.
library(emmeans)
emmeans(model2, ~ carsize*type) |>
contrast(interaction = "pairwise",
simple = "each",
combine = FALSE,
adjust = "bonferroni")
## $`simple contrasts for carsize`
## type = 1:
## carsize_pairwise estimate SE df t.ratio p.value
## 1 - 2 -20.000 4.67 30 -4.283 0.0005
## 1 - 3 50.833 4.67 30 10.886 <.0001
## 2 - 3 70.833 4.67 30 15.169 <.0001
##
## type = 2:
## carsize_pairwise estimate SE df t.ratio p.value
## 1 - 2 0.833 4.67 30 0.178 1.0000
## 1 - 3 52.500 4.67 30 11.243 <.0001
## 2 - 3 51.667 4.67 30 11.064 <.0001
##
## P value adjustment: bonferroni method for 3 tests
##
## $`simple contrasts for type`
## carsize = 1:
## type_pairwise estimate SE df t.ratio p.value
## 1 - 2 3.33 4.67 30 0.714 0.4808
##
## carsize = 2:
## type_pairwise estimate SE df t.ratio p.value
## 1 - 2 24.17 4.67 30 5.175 <.0001
##
## carsize = 3:
## type_pairwise estimate SE df t.ratio p.value
## 1 - 2 5.00 4.67 30 1.071 0.2928
A3_pg3
library(datarium)
se2 <- selfesteem2
se2_long <- se2 |>
pivot_longer(cols = c(t1, t2, t3),
names_to = "time",
values_to = "score")
se2_long <- se2_long |>
mutate(time = factor(time))
se2_long
## # A tibble: 72 × 4
## id treatment time score
## <fct> <fct> <fct> <dbl>
## 1 1 ctr t1 83
## 2 1 ctr t2 77
## 3 1 ctr t3 69
## 4 2 ctr t1 97
## 5 2 ctr t2 95
## 6 2 ctr t3 88
## 7 3 ctr t1 93
## 8 3 ctr t2 92
## 9 3 ctr t3 89
## 10 4 ctr t1 92
## # ℹ 62 more rows
Descriptives
describeBy(se2_long ~ treatment + time)
##
## Descriptive statistics by group
## treatment: ctr
## time: t1
## vars n mean sd median trimmed mad min max range skew kurtosis
## id 1 12 6.5 3.61 6.5 6.5 4.45 1 12 11 0.00 -1.50
## treatment 2 12 1.0 0.00 1.0 1.0 0.00 1 1 0 NaN NaN
## time 3 12 1.0 0.00 1.0 1.0 0.00 1 1 0 NaN NaN
## score 4 12 88.0 8.08 92.0 88.7 2.97 72 97 25 -0.75 -1.08
## se
## id 1.04
## treatment 0.00
## time 0.00
## score 2.33
## ------------------------------------------------------------
## treatment: Diet
## time: t1
## vars n mean sd median trimmed mad min max range skew kurtosis
## id 1 12 6.50 3.61 6.5 6.5 4.45 1 12 11 0.00 -1.50
## treatment 2 12 2.00 0.00 2.0 2.0 0.00 2 2 0 NaN NaN
## time 3 12 1.00 0.00 1.0 1.0 0.00 1 1 0 NaN NaN
## score 4 12 87.58 7.62 90.0 87.7 4.45 74 100 26 -0.41 -0.99
## se
## id 1.04
## treatment 0.00
## time 0.00
## score 2.20
## ------------------------------------------------------------
## treatment: ctr
## time: t2
## vars n mean sd median trimmed mad min max range skew kurtosis
## id 1 12 6.50 3.61 6.5 6.5 4.45 1 12 11 0.00 -1.50
## treatment 2 12 1.00 0.00 1.0 1.0 0.00 1 1 0 NaN NaN
## time 3 12 2.00 0.00 2.0 2.0 0.00 2 2 0 NaN NaN
## score 4 12 83.83 10.23 88.0 84.6 5.93 65 95 30 -0.62 -1.31
## se
## id 1.04
## treatment 0.00
## time 0.00
## score 2.95
## ------------------------------------------------------------
## treatment: Diet
## time: t2
## vars n mean sd median trimmed mad min max range skew kurtosis
## id 1 12 6.50 3.61 6.5 6.5 4.45 1 12 11 0.00 -1.50
## treatment 2 12 2.00 0.00 2.0 2.0 0.00 2 2 0 NaN NaN
## time 3 12 2.00 0.00 2.0 2.0 0.00 2 2 0 NaN NaN
## score 4 12 87.83 7.42 90.0 88.0 5.19 75 99 24 -0.47 -1.11
## se
## id 1.04
## treatment 0.00
## time 0.00
## score 2.14
## ------------------------------------------------------------
## treatment: ctr
## time: t3
## vars n mean sd median trimmed mad min max range skew kurtosis
## id 1 12 6.50 3.61 6.5 6.5 4.45 1 12 11 0.00 -1.50
## treatment 2 12 1.00 0.00 1.0 1.0 0.00 1 1 0 NaN NaN
## time 3 12 3.00 0.00 3.0 3.0 0.00 3 3 0 NaN NaN
## score 4 12 78.67 10.54 81.0 79.1 11.86 62 91 29 -0.39 -1.56
## se
## id 1.04
## treatment 0.00
## time 0.00
## score 3.04
## ------------------------------------------------------------
## treatment: Diet
## time: t3
## vars n mean sd median trimmed mad min max range skew kurtosis
## id 1 12 6.50 3.61 6.5 6.5 4.45 1 12 11 0.00 -1.50
## treatment 2 12 2.00 0.00 2.0 2.0 0.00 2 2 0 NaN NaN
## time 3 12 3.00 0.00 3.0 3.0 0.00 3 3 0 NaN NaN
## score 4 12 87.67 8.14 89.5 88.3 6.67 72 97 25 -0.69 -1.06
## se
## id 1.04
## treatment 0.00
## time 0.00
## score 2.35
Model
model <- aov_ez(id = "id",
data = se2_long,
within = c("treatment", "time"),
dv = "score",
type = 3,
anova_table = list(es = "pes"))
model
## Anova Table (Type 3 tests)
##
## Response: score
## Effect df MSE F pes p.value
## 1 treatment 1, 11 20.38 15.54 ** .586 .002
## 2 time 1.31, 14.37 7.24 27.37 *** .713 <.001
## 3 treatment:time 1.45, 15.90 6.06 30.42 *** .734 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
ANOVA
options(scipen = 999)
nice(model, es = "pes")
## Anova Table (Type 3 tests)
##
## Response: score
## Effect df MSE F pes p.value
## 1 treatment 1, 11 20.38 15.54 ** .586 .002
## 2 time 1.31, 14.37 7.24 27.37 *** .713 <.001
## 3 treatment:time 1.45, 15.90 6.06 30.42 *** .734 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG