20210717

library(tidyverse)

## -- Attaching packages -------------------------------------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## -- Conflicts ----------------------------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(ggpubr)
library(rstatix)

## 
## Attaching package: 'rstatix'

## The following object is masked from 'package:stats':
## 
##     filter

library(readxl)
options(scipen = 5)

1 基礎統計量と箱ひげ図

dat1 <- read_excel("C://Users/oldri/desktop/0718.xlsx", sheet = "Sheet2")
dat1$Group <- as.factor(dat1$Group)
dat1$Test <- as.factor(dat1$Test)
dat1$Group <- factor(dat1$Group, levels = c("Basic12","Basic34","Inter12","Inter34", "Inter56", "Advanced") )



dat1 %>% 
  group_by(Group, Test) %>%
  get_summary_stats(Score, type = "common")

dat_test1 <- filter(dat1, Test == "TEST1")
dat_test2 <- filter(dat1, Test == "TEST2")

ggboxplot(dat_test1, x = "Group", y = "Score", title = "TEST 1")

ggboxplot(dat_test2, x = "Group", y = "Score", title = "TEST 2")

2 1要因6水準

正規分布に従わない群があるため，ノンパラメトリックのKruskal-Wallis Testでやってみました。

2.1 TEST 1

2.1.1 外れ値

TEST 1 では外れ値がありません。

dat_test1 %>% 
  group_by(Group) %>%
  identify_outliers(Score)

2.1.2 正規性検定

Basic12, Inter12, Inter34 は正規分布に従わないので、Kruskal-Wallis Testを使います。

dat_test1 %>%
  group_by(Group) %>%
  shapiro_test(Score)

ggqqplot(dat_test1, "Score", facet.by = "Group")

dat_test1 %>% levene_test(Score ~ Group)

2.1.3 Kruskal-Wallis Test

群間の差が有意

res.kruskal <- dat_test1 %>% kruskal_test(Score ~ Group)
res.kruskal

2.1.4 Effict size

効果量が大

dat_test1 %>% kruskal_effsize(Score ~ Group)

2.1.5 Pairwise comparisons

#Dun法で多重比較する場合
pwc <- dat_test1 %>% 
  dunn_test(Score ~ Group, p.adjust.method = "bonferroni") 
pwc

#wilcos_testで多重比較する場合(オプション)
pwc2 <- dat_test1 %>% 
  wilcox_test(Score ~ Group, p.adjust.method = "bonferroni") 
pwc2

2.1.6 結果

pwc <- pwc %>% add_xy_position(x = "Group")
ggboxplot(dat_test1, x = "Group", y = "Score") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.kruskal, detailed = TRUE),
    caption = get_pwc_label(pwc),
    title = "TEST 1"
    )

2.1.7 まとめ

Dun法だと検出力がより厳しいようです。Pairwise comparisonsの節をご参照ください。

2.1.8 One-way ANOVA

おまけに一要因分散分析をやってみました。

res.aov <- dat_test1 %>% welch_anova_test(Score ~ Group)
res.aov

pwc <- dat_test1 %>% games_howell_test(Score ~ Group)
pwc

pwc <- pwc %>% add_xy_position(x = "Group", step.increase = 1)
ggboxplot(dat_test1, x = "Group", y = "Score") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
    )

# Welch One way ANOVA test
res.aov2 <- dat_test1 %>% welch_anova_test(Score ~ Group)
res.aov2

# Pairwise comparisons (Games-Howell)
pwc2 <- dat_test1 %>% games_howell_test(Score ~ Group)
pwc2

# Visualization: box plots with p-values
pwc2 <- pwc2 %>% add_xy_position(x = "Group", step.increase = 1)

ggboxplot(dat_test1, x = "Group", y = "Score") +
  stat_pvalue_manual(pwc2, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov2, detailed = TRUE),
    caption = get_pwc_label(pwc2)
    )

2.2 TEST 2

2.2.1 外れ値

dat_test2 %>% 
  group_by(Group) %>%
  identify_outliers(Score)

id 104 を除外しました。

dat_test2 <- dat_test2 %>% filter(!id == "104")

2.2.2 正規性検定

Basic34, Inter12, Inter34 は正規分布に従わないので、Kruskal-Wallis Testを使います。

dat_test2 %>%
  group_by(Group) %>%
  shapiro_test(Score)

ggqqplot(dat_test2, "Score", facet.by = "Group")

dat_test2 %>% levene_test(Score ~ Group)

2.2.3 Kruskal-Wallis Test

群間の差が有意

res.kruskal <- dat_test2 %>% kruskal_test(Score ~ Group)
res.kruskal

2.2.4 Effict size

効果量が大

dat_test2 %>% kruskal_effsize(Score ~ Group)

2.2.5 Pairwise comparisons

#Dun法で多重比較する場合
pwc <- dat_test2 %>% 
  dunn_test(Score ~ Group, p.adjust.method = "bonferroni") 
pwc

#wilcos_testで多重比較する場合(オプション)
pwc2 <- dat_test2 %>% 
  wilcox_test(Score ~ Group, p.adjust.method = "bonferroni") 
pwc2

2.2.6 結果

pwc <- pwc %>% add_xy_position(x = "Group")
ggboxplot(dat_test2, x = "Group", y = "Score") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.kruskal, detailed = TRUE),
    caption = get_pwc_label(pwc),
    title = "TEST 2"
    )

2.2.7 まとめ

Dun法だと検出力がより厳しいようです。Pairwise comparisonsの節をご参照ください。Test1と比べて，２つの多重比較の手法の違いがより少ないです。

2.2.8 One-way ANOVA

おまけに一要因分散分析をやってみました。

res.aov <- dat_test2 %>% welch_anova_test(Score ~ Group)
res.aov

pwc <- dat_test2 %>% games_howell_test(Score ~ Group)
pwc

pwc <- pwc %>% add_xy_position(x = "Group", step.increase = 1)
ggboxplot(dat_test2, x = "Group", y = "Score") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
    )

3 2要因

3.1 箱ひげ図

bxp <- ggboxplot(
  dat1, x = "Test", y = "Score",
  color = "Group", palette = "LocusZoom"
  )
bxp

3.2 外れ値

dat1 %>% 
  group_by(Test, Group) %>%
  identify_outliers(Score)

id 104 を除外しました。

dat1 <- dat1 %>% filter(!id == "104")

3.3 正規性検定

dat1 %>%
  group_by(Test, Group) %>%
  shapiro_test(Score)

ggqqplot(dat1, "Score", ggtheme = theme_bw()) +
  facet_grid(Test ~ Group)

box_m(dat1[, "Score", drop = FALSE], dat1$Group)

3.4 Two-way ANOVA

有意な交互作用が見られました。

res.aov <- anova_test(
  data = dat1, dv = Score, wid = id,
  between = Group, within = Test
  )
get_anova_table(res.aov)

習熟度の単純主効果

# Effect of group at each test
one.way <- dat1 %>%
  group_by(Test) %>%
  anova_test(dv = Score, wid = id, between = Group) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")

## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

one.way

# Pairwise comparisons between group levels
pwc <- dat1 %>%
  group_by(Test) %>%
  pairwise_t_test(Score ~ Group, p.adjust.method = "bonferroni")
pwc

Testの単純主効果

# Effect of time at each level of exercises group
one.way2 <- dat1 %>%
  group_by(Group) %>%
  anova_test(dv = Score, wid = id, within = Test) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")
one.way2

# Pairwise comparisons between time points at each group levels
# Paired t-test is used because we have repeated measures by time
pwc2 <- dat1 %>%
  group_by(Group) %>%
  pairwise_t_test(
    Score ~ Test, paired = TRUE, 
    p.adjust.method = "bonferroni"
    ) %>%
  select(-df, -statistic, -p) # Remove details
pwc2

# Visualization: boxplots with p-values
pwc <- pwc %>% add_xy_position(x = "Test")

bxp + 
  stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )

3.5 まとめ

二要因分散分析を行った結果，有意な交互作用が見られました。 習熟度の単純主効果を検定した結果，Test 1 では Basic 34 と Inter 12，Inter 12 と Inter 34， Inter 34 と Inter 56，Inter 34 と Advanced，Inter 56 と Advancedの差は有意ではありませんでした。Test 2では Basic 12 と Basic 34，Basic 34 と Inter 12，Inter 12 と Inter 34，Inter 56 と Advancedの差は有意ではありませんでした。